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ABSTRACT 

Pulsar observations provide a suite of tests to which stellar and binary evolutionary 
theory may compare. Importantly, the number of pulsar systems found from recent 
surveys has increased the statistical significance of pulsar population synthesis results. 
To take advantage of this we are in the process of developing a complete pulsar popu- 
lation synthesis code that accounts for isolated and binary pulsar evolution, Galactic 
spatial evolution and pulsar survey selection effects. In a recent paper we described 
the first component of this code and explored how uncertainties in the parameters of 
binary and pulsar evolution affected the appearance of the pulsar population in terms 
of magnetic field and spin period. We now describe the second component which 
focusses on following the orbits of the pulsars within the Galactic potential. In com- 
bination with the first component we produce synthetic populations of pulsars within 
our Galaxy and calculate the resulting scale heights as well as the radial and space 
velocity distributions of the pulsars. Correlations between the binary and kinematic 
evolution of pulsars are also examined. Results are presented for isolated pulsars, bi- 
nary pulsars and millisecond pulsars. We also test the robustness of the outcomes 
to variations in the assumed form of the Galactic potential, the birth distribution of 
binary positions, and the strength of the velocity kick given to neutron stars at birth. 
We find that isolated pulsars have a greater scale height than binary pulsars. This 
is also true when restricted to millisecond pulsars unless we allow for low-mass stars 
to be ablated by radiation from their pulsar companion in which case the isolated 
and binary scale heights are comparable. Double neutron stars are found to have a 
large variety of space velocities, in particular, some systems have speeds similar to 
the Sun. We look in detail at the predicted Galactic population of millisecond pulsars 
with black hole companions, including their formation pathways, and show where the 
short-period systems reside in the Galaxy. Some of our population predictions are 
compared in a limited way to observations but the full potential of this aspect will be 
realised in the near future when we complete our population synthesis code with the 
selection effects component. 

Key words: binaries: close - stars: evolution - stars: pulsar - stars: neutron - Galaxy: 
stellar content - Galaxy: kinematics and dynamics 



1 INTRODUCTION 

Detailed examinations of the stellar populations within 
the Galaxy, and indeed the greater Universe, have un- 
earthed many fascinating objects. Examples of such systems 
are gamma-ray bursts (Klebesadel, Strong & Olson 1973; 
Paczynski 1986; Bogomazov, Lipunov & Tutukov 2008), co- 
alescing double neutron stars (Rantsiou, Kobayashi, Laguna 
& Rasio 2008), X-ray binaries (Schreier et al. 1972; Liu, 
van Paradijs & van den Heuvel 2007; Galloway et al. 2008), 
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microquasars (Margon et al. 1979; Abell & Margon 1979; 
Combi, Albacete-Colombo & Marti 2008) and millisecond 
pulsars (Backer et al. 1982; Manchester et al. 2005). Cur- 
rent belief suggests that all the aforementioned systems arise 
from binary stars in which both stars are close enough to ex- 
perience a strong gravitational interaction with their com- 
panion. This may lead to mass transfer and depending upon 
the the mass ratio, types and ages of the two stars a plethora 
of different stellar and binary evolutionary phases may oc- 
cur. Recent surveys spanning a large range of observed wave- 
lengths now make it possible to monitor and analyse the 
combined properties of this variety of stellar populations. 
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This in turn can place further constraints upon our theo- 
retical understanding of stellar and binary evolution. For 
example, with the quickly increasing number of known low- 
mass stellar X-ray binary systems (187: Liu, van Paradijs 
& van den Heuvel 2007), high-mass X-ray binary systems 
(114; Liu, van Paradijs & van den Heuvel 2006) and pulsars 
(> 160(0; Manchester et al. 2005) it is possible to constrain 
features of neutron star (NS) and black hole (BH) forma- 
tion from their relative Galactic scale heights. These have 
been examined for NS and BH low-mass X-ray binaries (via 
observations, Jonker & Nelemans 2004) showing evidence 
that (contrary to previous belief) BHs may receive momen- 
tum from their formation mechanism - the supernova (SN) 
event. 

In this work our focus is on the Galactic population 
of pulsars. Observations of such objects occur primarily at 
radio (centimetre) wavelengths. These systems have intrin- 
sically weak signals (typically measured in mjy: Lorimer 
2005) and are observed as a regular series of pulses in time. 
Due to the frequency of the pulses we know these systems 
are compact (Hewish et al. 1968) , while the period derivative 
allows evolutionary models of these systems to be developed 
(e.g. Goldreich & Julian 1969; Ostriker & Gunn 1969; Gunn 
& Ostriker 1970; Bisnovatyi-Kogan & Komberg 1974; van 
den Heuvel 1984; Kulkarni & Narayan 1988; Chen & Rud- 
erman 1993) and tested (Dewey & Cordes 1987; Tauris & 
Bailes 1996; Dewi, Podsiadlowski & Pols 2005; Kiel, Hurley, 
Bailes & Murray 2008, hereafter Paper I). Finally, distances 
may be estimated from the dispersion of the pulse due to 
the electron density distribution within the Galaxy (Cordes 
& Lazio 2002) . From these observations we now believe that 
pulsars are magnetic rotating neutron stars and that the 
radio signal arises from the magnetosphere and is well colli- 
mated. Timing of these radio pulses gives the spin period P 
and the spin period derivative P from which the magnetic 
field and characteristic age of the pulsar may be inferred 
(see Paper I, and references therein for further details). 

Knowledge of pulsar space velocities provides con- 
straints on the effects of SNe on the NSs they give rise to. 
This may, in-turn, place constraints on the SN mechanism 
and NS structure. The birth velocities of pulsars (Lyne & 
Lorimer, 1994) combined with the work of Gott, Gunn & 
Ostriker (1970), Cordes, Romani & Ludgren (1993), Dewey 
& Cordes (1987) and Bailes (1989) demonstrated the neces- 
sity of asymmetric SN velocity kicks imparted on NSs by 
observing large isolated pulsar space velocities of order 1000 
km s~^. Later studies also found similar conclusions, again 
because observations suggest pulsar space velocities in ex- 
cess of the mean for normal field stars (Fryer & Kalogera, 
2001). Evidence for asymmetric SNe is also provided by the 
misalignment of binary pulsar spin vectors to the orbital an- 
gular momentum vector (e.g. Kaspi et al. 1996). The vectors 
would normally be expected to be coupled before the vio- 
lent SNe because of tidal effects and any occurrence of mass 
transfer (Bhattacharya & van den Huevel 1991; Hurley, Tout 
& Pols 2002). There is observational evidence of this mis- 
alignment for a number of pulsar binary systems including 
double NS systems (e.g. Kramer 1998). 

Another method in which we are able to constrain the 
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effects of SNe and to test the validity of our evolutionary as- 
sumptions (including the average magnitude of any kick de- 
livered by a SN) is to compare, with observations, the kine- 
matics of model pulsar populations within a model Galaxy 
(previous works include Dewey & Cordes 1987, Bailes 1989, 
Lorimer et al. 1993, Lyne et al. 1998 and Sun & Han 2004). 
There are two obvious pulsar populations we may recognise, 
those that are isolated pulsars and those pulsars within bi- 
nary systems. We may make further distinction with regards 
to the spin of a pulsar - those of a 'standard' spin period 
(P > Is), those of millisecond spin periods (P < 0.1 s: 
known as millisecond pulsars, MSPs) while those between 
these two period ranges are known as partially spun-up 
or partially recycled pulsars. The final pulsar population 
which we wish to point out here are those pulsars which 
reside in double compact binaries, in particular double NS 
binaries (such as PSR 1913 -I- 16: Hulse & Taylor 1975) 
and, within that population, double pulsar systems (PSR 
J0737 - 3039A&B: Burgay et al. 2003; Lyne et al. 2004). 

From a simple evolutionary analysis of these systems 
(see Paper I) you may expect there to be distinct scale 
heights for each pulsar population within the Galaxy. This 
is assuming that all pulsars are given an asymmetric SN 
kick velocity drawn from the same distribution (i.e. ignor- 
ing the possible electron capture SN scenario which may 
confuse this issue: see Paper I; Podsaidlowski et al. 2004; 
Ivanova et al. 2008). Isolated pulsars - it is reasonable to 
presume - would have a greater scale height than those pul- 
sars within binary systems. This is because binaries are on 
average heavier than single stars and also the binary orbit is 
able to absorb energy from the kick (as detailed in Hills 1983 
and Tauris & Takens 1998). Having been previously ejected 
from a binary or having always been a single star would 
allow the kick velocity to have greater effect on the stellar 
space velocity. Isolated pulsars that have felt two SN kicks 
during their lifetime (one indirect and one direct) would be 
expected to have the largest scale height of any pulsar pop- 
ulation (Bailes 1989). Say, for instance, that the first of two 
massive stars disrupts a binary system, giving momentum 
to both stars out of the Galactic plane, then the second 
massive star (now isolated) undergoes a SN and receives 
further momentum. Pulsars of different spin types, for ex- 
ample MSPs, can also be expected to have differing Galactic 
scale heights dependent upon the nature of their evolution- 
ary histories. If one assumes a MSP forms via Roche-lobe 
overflow mass-transfer from a low mass companion the sys- 
tem only passes through one SN event and this will occur in 
a close binary which can overcome larger velocity kicks to 
stay bound than, say, a standard pulsar in a wider (pre-SN) 
orbit (e.g. Bailes 1989; Portegies Zwart & Yungleson 1998). 
If, however, some MSPs are formed via wind accretion from 
a high-mass companion (as discussed in Paper I), which may 
itself form a NS, a fraction of the MSP population may feel 
two kicks. This may lead to an increase of the MSP Galac- 
tic scale height and also produce a method for isolated MSP 
production (as described by Narayan, Piran & Shemi 1991 
and Paper I). Finally, you may expect double pulsar sys- 
tems to have a greater Galactic scale height than the single 
pulsar binary systems because they feel two SNe. However, 
the fact that the binary system had to survive these two 
kicks requires the kick velocity imparted from both events 
to be small enough (or well directed) for the binary to sur- 
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vive. Therefore, although two kicks occur it seems plausible 
that double pulsar systems may not attain a greater scale 
height than their isolated cousins (Pfahl, Podsiadlowski & 
Rappaport 2005; Dewi, Podsiadlowski & Pols 2005). Fur- 
ther to this simple analysis, the relative ages of each system 
type will also play an important role in the observable scale 
height, because older systems will have had more time to re- 
lax (outwards) in the Galactic gravitational potential - the 
distributions diffuse over time. For a more in-depth review 
of pulsar evolution and kinematics see the living review of 
Lorimer (2008). 

The work described in this paper allows us to address 
these issues and to predict and compare the Galactic spatial 
distributions of pulsar populations. We follow on from our 
pulsar population synthesis in Paper I to not only evolve 
the stellar evolution of pulsars but move them within the 
Galactic gravitational potential. In other words we follow 
Galactic stellar, binary and kinematic evolution. As intro- 
duced in Paper I we are developing a code comprised of 
three modules; binpop (binary evolution), binkin (Galac- 
tic kinematics) and binsfx (synthetic survey simulations). 
An upcoming paper will describe the third module, binsfx, 
where we impose selection effects on the simulations, thus 
giving simulated data that can be compared directly to ob- 
servations. 

Section 2 outlines our binary evolution code (binpop) 
and details a necessary update to follow the evolution of a 
system that is disrupted owing to an asymmetric SN veloc- 
ity kick. In Section Owe describe the Galaxy kinematic code 
(binkin) which integrates the positions of pulsars - both in 
binary systems and isolated - forward in time within the 
Galaxy. This includes details of how the initial conditions 
for the Galactic population are chosen and the parameters 
involved. Results are given in Section |4] where, assuming a 
favoured binary and stellar evolutionary model of Paper I, 
we examine the pulsar population scale heights and veloc- 
ities that arise from different binkin model assumptions. 
This includes a detailed examination of the MSP-BH bi- 
nary population. In particular, we explore the formation and 
evolution of these systems. This is followed by a discussion 
of our findings and the main uncertainties involved in Sec- 
tion [S] 



2 RAPID BINARY EVOLUTION 

The first module, binpop, was described in detail in Pa- 
per I. Below we give an overview and also address the nec- 
essary modifications to binpop in order to correctly follow 
the Galactic positions of both members of a disrupted binary 
system. 

2.1 BINPOP 

binpop is a stellar/binary population synthesis package 
which convolves the binary stellar evolution (bse) code of 
Hurley, Tout & Pols (2002: hereafter HTP02) with realistic 
initial stellar and binary parameter distributions (as devel- 
oped in Kiel & Hurley 2006 and Paper I). Stellar evolution 
is included according to the formulae presented in Hurley, 
Pols & Tout (2000). Meanwhile, BSE attempts to account for 
all important binary evolutionary processes. These include 



tidal evolution, mass transfer, common envelope (CE) evolu- 
tion, stellar mergers, magnetic braking, orbital gravitational 
radiation and supernovae velocity kicks. In Paper I exten- 
sive additions were made to BSE, in terms of NS physics, so 
that pulsar evolution could be followed in detail. This means 
that aspects such as magnetic field decay, accretion induced 
field decay and spin-up, propeller evolution and pulsar death 
lines are now included. Inherent uncertainties in the vari- 
ety of binary and pulsar evolutionary processes requires a 
host of parameterised prescriptions to be incorporated into 
BSE. For example our lack of understanding of CE evolution 
is expressed as a parameter, acE, often referred to as the 
efficiency parameter. In terms of pulsar evolution there is 
the magnetic field decay time-scale, tb and the accretion in- 
duced decay time-scale, k, for example, which are uncertain. 
Over time the uncertainty in many parameters has decreased 
- albeit slightly - due to a large array of simulations and in- 
creasingly detailed observations (e.g. Lyne et al. 1998; Porte- 
gies Zwart & Yungelson 1998; HTP02; Belczynski, Kalogera 
& Bulik 2002; Podsiadlowski, Rappaport & Han 2003; Sun 
& Han 2004; Yusifov & Kucuk 2004; Hobbs, Lorimer, Lyne 
& Kramer 2005; Kiel & Hurley 2006; Cordes et al. 2006; 
Lorimer et al. 2006; Ferrario & Wickramasinghe 2007; Liu, 
van Paradijs & van den Heuvel 2007; O'Shaughnessy, Kim, 
Kalogera & Belczynski 2008; Belczynski et al. 2008 and Pa- 
per I). 



2.2 Binary evolution SN kick update 

When following the kinematic evolution of a binary sys- 
tem within the Galaxy we require knowledge of the Galactic 
gravitational potential - the acceleration felt on the binary 
centre of mass (CofM) owing to the Galaxy - as well as any 
internal sources of momentum that arise. The primary stel- 
lar evolutionary phase that can perturb an orbit or disrupt 
a binary system is a SN. If the SN occurs in a binary and 
enough material is ejected from the system during the event 
(more than half of the total binary mass) the binary may 
disrupt (Hills 1983). Along with the assumed instantaneous 
mass loss, if there is any asymmetry in the explosion (which 
is arguable in the case of BHs, see Podsiadlowski, Rappa- 
port & Han 2003 but note Pfahl, Podsiadlowski & Rap- 
paport 2005), the newly formed compact star will receive 
a velocity kick which the binary CofM will feel (Shklovskii 
1970; Lyne & Lormier 1994; Tauris & Takens 1998; HTP02). 
Depending upon the velocity kick direction and magnitude 
this may disrupt a binary or save it from mass loss-induced 
dissipation (Hills 1983; Kalogera 1998; Pfahl, Rappaport & 
Podsaidlowski 2003). 

The algorithm described by HTP02 allows realistic or- 
bital evolution modelling if the binary system survives the 
blast and stays gravitationally bound (eccentricity, e < 1). 
However, because binary systems cease to exist once they 
become unbound, HTP02 were not troubled with calculat- 
ing the recoil escape velocities of the two disassociated stars 
traveling on hyperbolic orbits with e > 1. Now that the 
Galactic spatial kinematics of both binary systems and iso- 
lated stars is a concern, knowledge of all associated velocity 
changes are required. To this end we formulate a disruption 
model (in Section [2.2.1|> and also describe similar methods 
derived by other groups (in Section r2.2.2|) . All methods con- 
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Figure 1. HTP02 orbital geometry as asymmetric SN occurs. 
Taken from Figure Al of HTP02. 




Figure 2. Orbital geometry of our disruption method after an 
asymmetric SNe. 



sidered here are generalised to allow for initially eccentric 
systems. 



2.2.1 BSE disruption model 

Before one can calculate the final velocities of the disrupted 
stars the binary system must be known to disrupt. We start 
by considering what is already within the capabilities of 
the rapid binary evolution code (cf. Appendix A of HTP02) 
which we outline here. This assumes a reference frame in- 
which the pre-SN CofM is at rest, Vs = 0, and the secondary 
star (the star not exploding) is at the origin (as shown in 
Figure The magnitude of the pre-SN relative orbital ve- 
locity is 



Vorb 



and vectorily is 

V — — Vorb (sin /3x + cos /3y) . 



(1) 



(2) 



The separation of the two stars is r = r[0, 1, 0] and (3 is 
the angle between r and V. Also, ^ = GMb where G is the 
gravitational constant, Afb = Mi + M2 and the subscripts 
denote the particular star (the primary star being 1 and the 



secondary 2). In this reference frame the two pre-SN stellar 
velocities are. 



Mb 



and 
V2 = 



Mb 



(3) 



(4) 



Furthermore, the orbital angular momentum J is expressed 



M1M2 
Mb 



r X V : 



M1M2 
Mb 



M1M2 



(5) 



where h is the specific angular momentum (aligned with 
the z-axis) and I = a{l — e^) is the semi-latus rectum. The 
primary star, the SN progenitor, is about to explode and 
receive a momentum impulse arising from a velocity kick. 



Vkick = Vkick (cos ijj cos (f)^ 4- sin u) cos (f)y -f sin 



(6) 



(where x, y and z are the unit directions vectors). The kick 
speed, Vkick, is modelled by a Maxwellian distribution. 



P (Vkic 



■ exp 



(7) 



as given by Hansen & Phinney (1997) with a dispersion V^. 
To facilitate any possible comparison to previous works, such 
as Hansen & Phinney (1997), Portegies Zwart & Yungelson 
(1998), HTP02 or Paper I, we take Kr = 190 km s"^ How- 
ever, a value of 265 km s"'^ has more recently been suggested 
by Hobbs et al. (2005) and there have also been suggestions 
of a bimodal kick distribution (Arzoumanian, Chernoff & 
Cordes 2002). The direction of the kick is specified by chos- 
ing two angles u) and 4> within the ranges < tj < 27r and 
— 7r/2 < (j) < -K 12 (as shown in Figure[T]|. 

Immediately after the asymmetric SN event the new 
velocity of the proto-NS is 



Vi = Vi + Vkic 



(8) 



and taking into account the instantaneous SN mass loss AM 
from the system we have a new relative velocity. 



V + Vk 



(9) 



The velocity of the new centre of mass relative to the old 
centre of mass is 



V, = 



Mns,, , AMMa,^ 

■ V kick + _ _ _ V 



Mi 



M'Mi, 



(10) 



(as in HTP02 equation A14). Following HTP02 the new ec- 
centricity, e , and semi-major axis, a , of the system can 
then be calculated. If this gives an eccentricity greater than 
unity then the system is disrupted and we need expressions 
for the runaway velocities of the stars. For further details on 
the SN treatment and the evolution of binary systems which 
survive the event see HTP02. 

Prior to the SN our chosen coordinate system had the 
orbital angular momentum vector directed along the z-axis 
but this will no longer be true for the post-SN system (unless 
Vkick ~ 0). However to simplify our post-SN calculations 
it is desirable to realign the vector with the z-axis. This 
realignment in the x-z plane, owing to the SN, is performed 
by a rotation, Rxz, around the y-axis such that the post-SN 
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unaligned orbital specific angular momentum, h = r x V = 
[0, r, 0] X [V"x, Vy, 14] becomes h" = [0, 0, K]. Note 
that we are using to denote the frame immediately after 
the SN and when referring to the coordinate system of 
the frame after rotation. Here the rotation is guided by u, 
the angle between the pre- and post-SN angular momentum 
vectors (see Equation A13 of HTP02). This rotation matrix 
allows us to map our final velocities back onto the original 
coordinate system. 

Now we consider the post-SN motion of the two stars 
which is governed by a hyperbolic conic section. We can 
no longer assume that the separation vector between the 
stars is aligned along the y-axis because the system is not 
necessarily at periastron. To account for this possible shift 
in coordinate system around the z-axis we calculate where 
in the orbit each star resides. Here we have 



- ' - ' e' \GM'r 



(11) 



with ip defined to be positive in the positive y-direction and 
zero along the positive y-axis. There are two possible hy- 
perbolic orbits for which each star may travel along - one 
in the positive y-region, the other in the negative y-region - 
which is governed by the direction of the y-component of the 
new relative velocity, that is, r ■ V . The sign of i/) depends 
upon the the sign of the r ■ V value, where i/; < when 
r V > 0. The post-SN binary mass must also be updated: 
Mb = AI-MS + M2 with Mns being the primary star mass. Us- 
ing this and the new semi-major axis we may now calculate 
the final velocities of the two stars in the pre-SN centre of 
mass reference frame. Assuming a velocity at infinity, Voo, 
which is directed along an asymptote of the hyperbolic orbit 
we have final velocities for the two stars of 



' _i M2 
Vlf = Vg — i?xz -tttVoo 

and 

A simple calculation gives the magnitude, 



Voo = 



gm; 



(12) 



(13) 



(14) 



The angle of the hyperbolic asymptote may be calculated 
from the angle a (as shown from Figure [2] which describes the 
post-SN coordinate system) where cos a = 1/e (as in Tauris 
& Takens 1998) and a is always positive. This restricts cr to 
range from — > n/2. With our two angles we define the 
difference angle y = a — which is used in rotating the 
coordinate system around the z-axis. Separating Vif and 
V2f into component form gives us: 



r Vkick cos LO COS (f> + ^ ^, ^ ^ Vorb Sm fj 

— Voo COS u sin 7, 



V 



— Vkick sm UJ COS (p + ^ , , Vorh COS [i 



-Voo COS 7, 



Vu 



Mns 
Mi 



Vkick sin ( 



(15) 



(16) 



(17) 



V2i: 



V2ty = 



and 



V2{: 



Mns,, , , AMM2.. . ^ 

— Vkick COS LO COS (f> H ^ ^, ^ ^ Vorb Sm /£> 



m; 

+V00 COS f sin 7, 



M;Mb 



Mns,, . ^ , AMM2 ,, 

— Vkick sm LJ COS (f> -\ ; Vorb COS p 



K 

+V00 cos 7 



Mns , , 

—-^ Vkick smt 

Mk 



(18) 



(19) 



(20) 



We also account for coalescence of the two stars if the 
newly formed compact star velocity kick is directed towards 
the companion. Coalescence occurs if the companion radius, 
R2, is greater than periastron or the distance of closest ap- 
proach: if i?2 > a (e — 1) we assume the stars coalesce (see 
also Tauris & Takens 1998; Belczynski et al. 2008) and the 
merger outcome has the final centre of mass velocity, 
and the mass of the NS and companion are combined (see 
HTPQ2 for further details of merger outcomes). 

Equations 1151 to 1201 are the velocities calculated within 
the BSE kick routine and are communicated into binkin. 
Before adding these to the Galactic binary centre of mass 
velocity at the time of SN it is first necessary to perform 
a random orientation of the pre-SN binary orbital plane, 
which until now has been fixed in the Galactic xy-plane. We 
randomly choose three Euler angles ob, /3b and 7e, within 
the ranges ^ aE,7B < 27r and ^5 /3b < tt, to give a 
3D rotation of the velocities in Equations [T^] and 1131 The 
post-SN velocities (pre-SN CofM velocity plus the rotated 
disruption velocities) are then used within the kinematic 
routine to calculate the subsequent velocities and positions 
within the Galaxy of the pulsar and its former companion. 

2.2.2 Related disruption models 

There are two other groups who have independently devel- 
oped models for deriving the run-away velocities of stars 
from a disrupted binary. The method published recently by 
Belczynski et al. (2008) is similar to our demonstration and 
is also generalised for arbitrary eccentricity. The main dif- 
ference is that it allows a velocity Vimp to be imparted to the 
secondary from the expanding shell of the SN. In our work 
we essentially assume that Vimp ~ which has minimal ef- 
fect except in some cases of small pre-SN orbital separation 
(Kalogera 1996; Belczysnki et al. 2008). Tauris & Takens 
(1998; hereafter T&T98) have developed a relatively sophis- 
ticated model. The major difference between our model and 
the T&T98 model is the coordinate system scheme and the 
latters assumption that the companion star may have mo- 
mentum imparted directly onto it from the SN blast wave 
(as in Belczynski et al 2008 and Dewey & Cordes 1987). The 
companion star may also have some fraction of mass stripped 
off it and/or ablated owing to the impact of the shell of ma- 
terial ejected from the primary. To include the possibility 
of investigating the effect of these additional considerations 
we have worked through the T&T98 demonstration and im- 
plemented it as an option in BSE, generalised to eccentric 
orbits. However, we do not exercise this option in this work. 
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Figure 3. The speeds of the two stars following the SNe. The 
thick line represents the NS recoil speed distribution while the 
thin line is the companion recoil speed distribution. Included is 
the assumed asymmetric SN kick distribution assuming a disper- 
sion of 190 km (dashed line). 



2.2.3 Disruption model illustration 

To illustrate the effect the binary orbit has on the runaway 
velocities of disrupted stars we produce a simple population 
of binary systems. For this population the primary mass, 
Ml, is randomly selected from a flat distribution ranging 
from 10 — 20 Mq, the secondary mass, M2, from a flat dis- 
tribution ranging from 0.1 — 20 Mq and the orbital sepa- 
ration is selected randomly from a flat distribution from 1 
to 10 000 -R0. All systems are initially circular to simplify 
the analysis. The radius of the secondary star is linked to 
its mass by R2 = 1.3M2 '® if M2 is greater than unity and 
by i?2 = M2 otherwise. We make sure the system is not 
in contact at birth. We then let the primary undergo a SN 
that leaves a NS with Mns = 1-4 Mq and assume the rem- 
nant is given a kick from Equation [7] with Va = 190 km s~^ 
(in accordance with Hansen & Phinney 1997). The post-SN 
velocities for disrupted stars are calculated using the BSE 
method detailed above. 

After generating a million systems we flnd that the ma- 
jority (99%) become unbound and the incidence of coales- 
cence is negligible. In Figure [3] we see the distributions of 
NS and companion star recoil speed and compare this to 
the Vkicii distribution, i.e. the distribution for a popula- 
tion of standard single NSs. The first item we wish to note 
is the difference in the typical velocities received by both 
stars: the NS, which directly experiences the additional mo- 
mentum imparted from the asymmetric SN, will most likely 
depart the binary system with a greater velocity than the 
companion star (relative to the CofM). The second point of 
interest is the similarity between the NS recoil speed dis- 
tribution and the kick distribution. Clearly not all of the 
momentum imparted onto the NS goes into the NS recoil 



velocity, some of the momentum is instead transported into 
the CofM momentum, consumed by the disruption of the 
binary system and converted into additional velocity of the 
companion star. Therefore, although the shape of the NS 
recoil speed distribution is consistent with the kick distribu- 
tion the NS distribution is shifted somewhat to lower values. 
In regards to this, observational pulsar velocity studies that 
do not account for the possibility of a fraction of the sam- 
ple being disrupted from binary systems may underestimate 
the underlying SN kick distribution. It also suggests a pos- 
sible mechanism for any bimodality found in the velocity 
structure of pulsar observations - similar to that found by 
Arzoumanian, Chernoff & Cordes (2002) from which they 
concluded a bimodal asymmetric SN kick distribution, or 
possibly binary disruption effects, could cause such detected 
velocities. We note that the form of the underlying orbital 
period (or separation) distribution of the model binaries will 
affect the distribution of recoil speeds - with an increased 
proportion of short-period systems leading to an increased 
difference between the NS recoil speed and kick distributions 
- and we have not explored this aspect in detail here. 



3 GALACTIC KINEMATICS 

3.1 Galactic gravitational potentials 

Much work over the years has lead to estimates of the Galac- 
tic gravitational potential. Miyamoto & Nagai (1975) gen- 
eralised the work of Toomre (1963) who calculated flattened 
Plummer (1911) models for the Galaxy. Since then further 
observations have lead to estimates by Carlberg & Innanen 
(1987) of disk-halo, bulge and nucleus potentials which in 
turn have been updated by Kuijken & Gilmore (1989). Kui- 
jken & Gilmore (1989) used more extensive observations of 
Galactic stellar densities, which allows the mapping to an 
assumed (to first order) smooth time-independent galactic 
gravitational potential. 

The KG89 model potential is, 



where. 



(21) 
(22) 



=(23) 



(24) 
(25) 



and the parameter values for each region are: 



disk/halo: /3i = 0.4, 132 = 0.5, /Ja = 0.1, hi = 0.325 kpc, 
h2 = 0.090 kpc, h-i = 0.125 kpc, a = 2.4 kpc, b = 5.5 kpc, 
Mdisii = 1.45 X 10" Me 
nucleus: b = 0.25 kpc, Mnuc = 9.3 x lO'' Mq 
bulge: b = 1.5 kpc, Mbuigc = 1-0 x 10^° Mq. 

If necessary, use of the KG89 model allows us to com- 
pare results to previous binary pulsar population synthesis 
works such as Lorimer et al. (1993). While now considered 
somewhat outdated, the KG89 potential is still in use within 
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recent observational works, such as, Freire, Ransom & Gupta 
(2007) who use it in their calculation of the pulsar spin pe- 
riod when accounting for observational effects of the acceler- 
ation between the Solar System Barycentre and NGC 1851 
- the globular cluster in which their observed pulsar resides. 
However, a recent appraisal of the most promising Galactic 
gravitational potentials to use was completed by Sun & Han 
(2004). Their favoured method for ease of implementation 
is that given by Paczynski (1990: see below). Sun & Han 
(2004) also commented favourably on the work of Dehnen 
& Binney (1998) who fit a multi-parameter mass model to 
kinematic data of the Milky Way. Sun & Han (2004) found 
that the Dehnen & Binney (1998) model is overly compli- 
cated to set up and manipulate, while the simpler Paczynski 
(1990, hereafter P90) model is as accurate as the Dehnen & 
Binney (1998) model. As such we also include the P90 model 
in our work. 

The P90 model, like the KG89 model, follows the po- 
tential of Miyamoto & Nagai (1975) for the disk (i = 1) 
and spheroid components {i — 2). Equation 9 of Paczynski 
(1990) is, 

GMi 



<E-f(i?,2) = 



rR2 



-t-62){l/2) 



|2\(l/2) ' 



(26) 



where R — y^x^ + y'^). The P90 model, however, differs 
from the KG89 model not only with the assumed constant 
values used within Equation [26] (for the disk potential) but 
also with the assumed form of the halo potential. 



^P, , GMh 



1 , r , rh f r 

- ln(l + —) H atan I — 



(27) 



where = + is used to simplify the equation. The 
parameters being: 

disk (i = 1): ai = kpc, bi = 0.277 kpc. Mi = 1.12 x 



*bulgo(»') 

and 



GM, 



bulge 



r + Co 



/ N GA/vir , 

*NFw(r) = — t;!—^^ 



\ / vir 



(29) 



(30) 



Here Gx = ln(l -I- Gx) — Gx/(1 + Gx), while the values used 
within Xue08 vary depending upon the assumed halo de- 
scription. We also note here that the form of $npw(?') is 
exactly that of Smith et al. (2007), who provide differing 
values for the virial mass Afvir, radius r^h and concentra- 
tion Gx. XueOS match their observed circular Galactic stellar 
velocity estimates to smooth particle hydrodynamical sim- 
ulations from which they provide values for Afvir, r^ir and 
Gx. The values from Xue08 used in our work are: 

disk: Afi^.k = 5 x 10^° Afg, b = 4 kpc, 
bulge: Af^^i_ = 1.5 x 10^° Mq,Co = 0.6 kpc. 



= 278 kpc, Gx = 11.. 



These Galactic gravitational potential models are now 
a part of binkin, updating the original algorithm based on 
Lorimer, Bailes & Harrison (1997), and extending upon sim- 
ilar population synthesis models such as that of Faucher- 
Giguere & Kaspi (2006). 

Figure |4] depicts the three assumed Galactic gravita- 
tional models. In particular we wish to point out the inner 
region of the Xue08 model, which contains a smaller restor- 
ing force than the other two models. We return to this later 
in Section m The KG89 gravitational potential model decays 
faster than the other two models beyond the central Galac- 
tic region and the P90 model rotational curve follows the 
KG89 model in the inner region of the Galaxy while beyond 
a Galactic radius of ~ 20 kpc the rotational curve flattens 
off similarly to the XueOS model. 



spheroid (i = 2): a2 = 3.7 kpc, 62 = 0.20 kpc, M2 = 
8.07 X 10^" A4q, 
halo (i = h): n, = 6.0 kpc, Afh = 5.0 x 10^" Mq. 

The KG89 and P90 models are both based on old obser- 
vations of the stellar neighbourhood. As such, these models 
are only considered accurate out to a radii of ~ 12 kpc. More 
recent observations completed in the Sloan Digital Sky Sur- 
vey (SDSS: York et al. 2000), within the Sloan Extension 
for Galactic Understanding and Exploration (SEGUE: Lee 
et al. 2008) program, have allowed the Galactic gravitational 
potential to be measured out to a radii of ~ 60 kpc (Xue 
et al. 2008). To do this Xue et al. (2008) have made line 
of sight velocity measurements of ~ 2500 blue horizontal 
branch stars which are converted into circular velocity es- 
timates of the Milky Way. Ultimately the work of Xue et 
al. (2008) is completed to probe the halo of our Galaxy and 
thus they do not examine in any detail the inner Galactic 
potential. However, at this stage we consider their complete 
assumed Galactic gravitational potential as an option in our 
work. The Xue et al. (2008, hereafter XueOS) model makes 
use of the following exponential disk, Hernquist (1990) bulge 
and Navarro, Frenk & White (1996; NFW) halo potentials 
respectively: 



GA/i^,; 



[1 



3.2 Initial conditions and integration method 

Once we have our assumed Galactic potential, which is in 
cylindrical coordinates, <1> (r, z), the progenitor pulsar sys- 
tems must be given some initial Galactic position, -Rinit, se- 
lected randomly from a given distribution. The most straight 
forward distribution is to assume a thin disk with some max- 
imum height and radius. This simple method allows the sys- 
tem to relax over time into a similar distribution to the 
observed stellar number density distribution in height with 
respect to the plane (see distributions given in Sun & Han 
2004). However, as shown in Section [4.11 this overestimates 
the number of systems (in this case pulsars) found in the 
central region - the deficiency of observed pulsars within 
the central region is believed to be a real lack of pulsars and 
not caused solely by observation selection effects (Lorimer 
et al. 2006). However, Lorimer et al. (2006) caution readers 
that more observations are required to provide a definitive 
result. 

To combat this over density within the central regions 
a preferred option is to use the distribution of SN remnants 
developed in Paczynski (1990): 

P{Rinit)dR = an (i?init/i?exp) exp (-i?init/i?cxp)d-R , (31) 



exp 



-r/bl 



(28) 



where Rc 



4.5 is simply an exponential scale length of 



the radial distribution and aa is a constant of integration 



8 



P.D. Kiel and J.R. Hurley 



^2 




equations of motion to evolve the system position forward 
in time. These are (Paczynski 1990): 
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Figure 4. The three Galactic gravitational models arc depicted 
in two different manners here to illustrate their properties and dif- 
ferences. The top panel (a) shows the circular velocity for a range 
of Galactocentric radial positions in the plane of the Galactic 
potential. The bottom panel (b) depicts modulus acceleration in 
increasing height from the plane. 



equal to 1.0683 over i? = ^ 20 kpc. A third option is 
to use the distribution derived by Yusifov & Kucuk (2004; 
their Equation 17) from observations of OB type Population 
I stars: 



P(i?init)d-R = an, (-Rinit/i?©) exp 



Rr.; 



dR. 



(32) 



If OB stars are assumed to be the progenitors of NSs then 
this can be taken as the Galactic pulsar progenitor birth ra- 
dial distribution. Here hr ~ /24: ~ 606 (for ^ Anit ^ 
20 kpc). We utilise all three radial distributions - thin disk, 
Paczynski (1990) and Yusifov & Kucuk (2004) - in our 
models to generate birth locations. However, because the 
Paczynski (1990) distribution has had much use in the past 
we assume this to be our standard description. 

In terms of the initial distribution of systems in height 
from the plane we simply assume a uniform distribution 
with maximum height of |2;maxi|. Systems over time relax 
outwards in \z\ and such a simple initial \z\ distribution 
compares well to those favoured in Sun & Han (2004). Fur- 
thermore, according to Paczynski (1990) as long as the sys- 
tems are born relatively close to the Galactic plane (few 
hundred parsec) the initial distribution in \z\ for energetic 
populations is redundant. In the future the initial spatial 
distribution will contain spiral arms, similar to that com- 
pleted within Faucher-Giguere & Kaspi (2006) who suggest 
that because Galactic arm structure is visible in large ob- 
servational surveys it is necessary for any realistic pulsar 
population synthesis simulation to model this structure. 

Once a position is found for a system an initial velocity 
is simply calculated from the (estimated) Galactic gravita- 
tional potential at that point. With knowledge of the posi- 
tion and space velocity the next step is to solve four coupled 



dR 
dt 
dK 

dt 



Vr, 



dz 
It 

dz 



dVR 
dt 



(33) 



These four equations are found by calculating the accelera- 
tion in R, <j) and z (Vr and 14 are velocities in R and z re- 
spectively) induced onto a test particle by the gravitational 
potential. Assuming an axisymmetric potential around z 
produces constanst angular momentum, L, felt by a test 
particle. To integrate forward in time a fourth order Runge- 
Kutta integration routine is used, similar to that used by 
Paczynski (1990) and Lorimer et al (1993). 

It is now possible for us to evolve the complete Galactic 
orbital evolution of a system of interest. If a SN event occurs 
and the system is not disrupted the velocity injected into the 
system by a SN, calculated within BSE, is simply vectorially 
added to the known Galactic velocity of the system centre 
of mass at the time of the SN. If the system is disrupted 
the run away velocities of the two stars - as calculated in 
Section [2.21 - are, again, vectorially added to their previous 
Galactic velocity (that of their system of origin) . In this way 
we are able to follow the complete Galactic orbital history of 
a system (star or binary), even if it passes through two SNe 
and with (or without) binary disruption. Because we assume 
no interaction between orbiting systems we are able to evolve 
each system separately, one after another (or in parallel). A 
beneficial consequence resulting from this assumption is that 
BINKIN is faster to run than other dynamical codes, such as 
typical N-body codes (McGlynn 1984). Of course, if one is 
dealing with compact stellar clusters dynamical interactions 
between the stellar components are very important to follow. 



4 PULSAR POPULATION STATISTICS 

We now describe the results of a series of population syn- 
thesis calculations that utilise our binpop and binkin mod- 
ules to follow the stellar/binary and kinematic evolution of 
a population of binary stars to produce artificial Galactic 
pulsar populations. The primary aim of this section is to 
predict scale heights and other kinematic characteristics for 
populations of pulsars, firstly assuming that all pulsars can 
be detected. We also compare the kinematics of our model 
pulsar populations with available kinematic tracers found 
from pulsar survey observations (such as those produced in 
Yusifov & Kucuk 2004; Lorimer et al. 2006 and Hobbs et al. 
2005). For now we hold back from making direct statistically 
significant comparisons as this requires modelling of selec- 
tion effects which will be covered in detail in a companion 
paper (the binsfx component as mentioned in Section 
Here we focus more on showing how modifying certain pa- 
rameters affects the final pulsar population kinematics, in 
terms of scale heights and space velocity distributions. 

The first step is to evolve a population of binaries 
within BINPOP. For this we proceed using our favoured model 
from Paper I (Model Fd). This sets choices for binpop stel- 
lar and binary evolutionary parameters of: solar metallic- 
ity Z = 0.02; a maximum possible NS mass of 3 Mq and 
acE = 3. It also sets the following parameters governing 
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Table 1. A summary of the BINKIN models used within this paper 
and their main assumptions. The first column provides a label for 
each model. This is followed by the assumed Galactic radial dis- 
tribution of pulsar birth locations, -Rinit, where Paczynski stands 
for the distribution suggested by Paczynski (1990; see our Equa- 
tion |3TJ, Flat is the flat distribution described in Section l3.2l and 
Yusifov & Kucuk is the birth distribution suggested in Yusifov & 
Kucuk (2004; see our Equation I32|l. The third column gives the 
Galactic gravitational model used (see Section [XTJ. The next col- 
umn is the maximum height with respect to the Galactic plane, 
l^maxl, that we consider for our pulsar distribution and scale- 
height calculations, followed by our assumed value of the SN kick 
dispersion, Va- Note that Model G evolves a population of single 
stars while all other models start with a population of binaries. 



Model 


-Rinit distribution 




-^max 1 




A 


Paczynski 


P90 


10 kpc 


190 


B 


Flat 


POO 


10 kpc 


190 


C 


Yusifov & Kucuk 


POO 


10 kpc 


190 


C2 


Yusifov & Kucuk 


POO 


2 kpc 


190 


C3 


Yusifov & Kucuk 


POO 


20 kpc 


190 


D 


Paczynski 


KG89 


10 kpc 


190 


E 


Paczynski 


XueOS 


10 kpc 


190 


F 


Paczynski 


POO 


10 kpc 


550 


G 


Paczynski 


POO 


10 kpc 


265 



pulsar evolution: re — 2000 Myr; k = 3000; no propeller 
evolution; the initial pulsar period and magnetic field pa- 
rameter selections linked to the strength of the SN velocity 
kick; the angular momentum accreted by the pulsar is vari- 
able; no electron capture SNe; and no beaming of pulsars 
(see Section [2. II and Paper I for details). SN kicks using the 
BSE prescription are given to NSs (see Section |2.2|) . while 
to keep the required number of models down to a minimum 
we only use the curvature radiation death line model (Hard- 
ing, Muslimov & Zhang 2002). Unless otherwise stated we 
take 10^ binary systems with initial parameters selected in 
the same manner as within Paper I and with the limits of 
5 — 80 Mq for primary mass, 0.1 — 80 Mq for secondary 
mass and 1 — 30 000 days for orbital period. The Galactic 
age is assumed to be 10 Gyr. Each binary is evolved to this 
age and for those that create pulsars the evolution history 
e.g. SN occurrence times and velocities, is saved as input for 

BINKIN. 

The next step is to take each of the binpop binaries and 
follow their corresponding kinematic evolution in binkin. 
For each binary a random birth age is assigned and the evo- 
lution followed from this age up to the age of the Galaxy. For 
this we begin by defining a standard model which we will 
call Model A. This uses the Paczynski (1990) distribution 
for setting the initial Galactic radial positions of the bina- 
ries (see Equation [31} with a maximum initial height off the 
plane of |zmaxi| ~ 75 pc. It also assumes the P90 form of 
the Galactic potential (see Equations [26] and I27|l and sets 
Vo- — 190 km s~^ (used within binpop) as the dispersion of 
the SN velocity kick distribution. Further models arise due 
to variations of these choices and are listed in Table [1] 

For each model we examine the scale heights for a range 
of pulsar systems. These are given in Table (2] We take the 
scale height to be that distance in |z| for which the number 
of stars within that distance is 63% (~ twice the e-folding 
distance) of the entire population. The most prolifically ob- 



served pulsar system is what is known as a standard pulsar. 
Here we define a standard pulsar as one which satisfies 

logS ^ -2.5 X logP + 8.1. (34) 

This equation artificially divides the 'standard' pulsar 'is- 
land' from all other radio pulsars in the B — P diagram 
(see Paper I). We define a MSP to be a pulsar spinning 
more rapidly than P — 0.02 s. All other pulsars bridge these 
two pulsar types - islands within the PP plane (see also 
the description given in Paper I). We also distinguish be- 
tween binary and isolated pulsars. It is possible to compare 
our model results in a limited manner to observations. To 
do this we make use of the ATNF Pulsar Catalogue which 
provides us with ~ 1610 Galactic plane pulsars (we ignore 
pulsars from the catalogue that have any association with 
another object, for example with a globular cluster or exter- 
nal galaxy). Approximately 1550 of these are isolated. Only 
15 standard pulsars are found in binary systems within the 
Galactic disk. Of the total observed pulsars there are 65 
MSPs of which 19 are isolated. We show the scale heights of 
the catalogue pulsars within Tabled 

The radial distributions of Galactic pulsars resulting 
from our set of models are shown in Figure [5] (left-hand side 
panels). We also compare a subset of the models in more 
detail in Figure [§] and include a comparison to the initial 
distributions used in the models and also the pulsar distri- 
bution suggested by Yusifov & Kucuk (2004) which is based 
on observations. We also explore the pulsar population 3D 
space velocity distributions. These are shown in the right- 
hand side panels of Figure [5] for the models in Table [1] 

Our results are analysed in more detail in four parts. 
In Section [4. II we examine the effect of varying the assumed 
pulsar birth radial distribution. This analysis makes use of 
Models A, B and C. Within Section 14.11 we also consider 
how modifying the target area considered (the 'observable' 
Galactic area) in our scale height calculations affects the 
scale height values of pulsars produced in Model C. This is 
completed with the use of Models C2 and C3. Section 14.21 
analyses different forms of the Galactic gravitational poten- 
tial by comparing Models A, D and E pulsar scale heights, 
final radial distributions and final velocity distributions. We 
then consider the effect of varying the assumed SN velocity 
kick distribution with Models A, F and G in Section 14.31 
Finally, after examining differences in bulk pulsar proper- 
ties, we explore in detail the MSP population of Model C 
in Section 14.41 We make use of our MSP analysis to further 
investigate the effect of model assumptions such as the ini- 
tial scale height. Galactic age and the number of systems 
evolved. 

4.1 Initial distributions and target area 

To begin with we focus on Model A. We note that unless oth- 
erwise specified, when calculating scale heights we consider 
only pulsars located within the Ga lactic region defined by 
1^1 ^ l^maxl kpc and r — y'x^ + y^ ^ 30 kpc. Firstly com- 
paring the scale heights of the isolated and binary pulsars 
we see that as expected (see Section [1} it is the isolated pul- 
sars which have the greatest scale height when considering 
all pulsars. It is also no surprise that the scale height of the 
complete population (top row of Table [2} is closely aligned 
with the isolated pulsar scale height. This is because 91% 
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Table 2. Model scale heights for a range of pulsar types in the Galaxy. 







Model A 


B 


C 


D 


E 


F 


C2 


C3 


Type 




Both 


1.39 


1.49 


1.25 


1.33 


2.58 


1.93 


0.58 


1.40 


All 


Isolated 


1.43 


1.53 


1.30 


1.37 


2.64 


1.93 


0.59 


1.46 




Binary 


0.96 


1.10 


0.79 


0.91 


1.98 


2.00 


0.46 


0.83 




Both 


1.43 


1.53 


1.30 


1.36 


2.63 


1.93 


0.59 


1.45 


Standard 


Isolated 


1.43 


1.53 


1.30 


1.37 


2.64 


1.93 


0.59 


1.46 




Binary 


0.73 


0.85 


0.59 


0.71 


1.68 


1.58 


0.38 


0.61 




Both 


1.00 


1.15 


0.82 


0.95 


2.04 


2.02 


0.47 


0.88 


All MSPs 


Isolated 


1.76 


1.67 


1.76 


1.49 


2.90 


2.33 


0.60 


2.21 




Binary 


0.99 


1.14 


0.82 


0.94 


2.03 


2.02 


0.47 


0.87 



of pulsars in Model A are isolated at the end of the simula- 
tion (even though all stars are initially in binaries). When 
considering only standard pulsars the domination of the iso- 
lated component is even greater: 99% of standard pulsars 
are isolated. However, the tables are turned when we look 
at the MSP population: the binary MSP population makes 
up 99% of all MSPs within Model A (a greater percentage 
than what is observed, however, see Section |4AT] for further 
discussion on this point). These relative numbers explain 
why in Table [2]the scale height of standard pulsars is insen- 
sitive to standard binary pulsars, and the same can be said 
for the total MSP scale height compared to isolated MSPs. 

For the binary pulsars we find that the standard pulsar 
population has a lower scale height than for the MSP pop- 
ulation. This suggests that on average binary MSPs receive 
greater post-SN recoil velocities than their standard binary 
pulsar counterparts. The reasons for this were alluded to in 
Section [1] but we reiterate them here (see also the findings 
of StoUman & van den Heuvel 1986; Bailes 1989). Basically, 
the very existence of an MSP relies on the occurrence of 
mass-transfer on to the NS which in turn requires a close 
binary. Such a binary will have a greater binding energy 
than the equivalent standard (wider) binary pulsar and can 
therefore survive a greater SN velocity kick as there is more 
energy to overcome for disruption. As a result it is possible 
for the proto-binary MSP system to be given a faster recoil 
velocity (with respect to the systems initial CofM). The iso- 
lated MSP population scale height also increases compared 
to the entire isolated pulsar population. Again this is not 
surprising given the model isolated MSP formation scenario 
as addressed in Paper I and within Section [T] Owing to the 
SN that produced the NS progenitor of the MSP, and the 
SN of the companion star that disrupted the system, iso- 
lated MSPs can receive greater recoil velocities than any 
other pulsar population - hence the largest population scale 
height (this formation mechanism is also discussed further 
in Section 1331 . 

To measure what effect the initial Galactic birth distri- 
bution has on the final model pulsar population we compare 
Models A, B and C. These models assume the birth distri- 
bution of Paczynski (1990: from the observed distribution of 
SN remnants), a uniform thin disk, and the birth distribu- 
tion of Yusifov & Kucuk (2004: derived from observations of 
OB stars), respectively. The final radial distributions (at a 
Galactic age of 10 Gyr) for the pulsars in these models are 



compared in the upper-left panel of Figure[S] Both Models A 
and C have distributions that peak away from the Galactic 
centre (refiecting their initial distributions). However, the 
peak of Model A is shallower and the distribution is more 
extended than Model C. We note that the number of pulsar 
systems 'observed' in Model C is higher than in Model A 
(by a factor of ~ 1.1). For Model B we find that the pul- 
sar radial distribution peaks at the Galactic centre and then 
decays with increasing radius. The space velocity distribu- 
tions for the pulsars in the three models are compared in 
the upper-right panel of Figure [5] and we see that there is 
no discernible difference. 

The relation of the final pulsar radial distribution to 
the assumed birth distribution of binaries can be seen in 
Figure |6] for Models A and C. Here all distributions are nor- 
malised to unity to aid comparison of the peak position and 
distribution shapes. We see that the birth distribution of 
Model A is broader than for Model C and this is refiected 
in their final shapes. However, the width of the distribution 
increases with time in both cases, while the peak of the dis- 
tribution moves in towards the Galactic centre, which is a 
typical effect of Galactic potentials (Sun & Han 2004). Ini- 
tially Model A peaks at a radius of 4.5 kpc which moves 
inwards to a radius of 3.9 kpc at 10 Gyr. For Model C the 
peak moves from 5.0 kpc to 4.5 kpc. 

In Figure |S] we also compare the model distributions to 
the distribution of an observed sample of pulsars presented 
by Yusifov & Kucuk (2004). We note that at this stage we 
are not including selection effects in our models so a direct 
comparison with observations is not possible. However, com- 
parison with the Yusifov & Kucuk (2004) sample, which in- 
cludes selection effects somewhat by being limited to pulsars 
with P > lO"^'^ s s~^, can still provide a meaningful guide 
to discerning between our models. Although not included in 
Figure |6] we can see immediately that Model B is an unreal- 
istic model of the Galactic pulsar population. The apparent 
deficit of observed pulsars in the inner region of the Galaxy 
can not be reproduced by assuming all binaries are born in 
a uniform thin disk - we require a paucity of pulsars to be 
born in the central region of the Galaxy when attempting to 
match observations (similar to Paczynski 1990; Sun & Han 
2004). This lack of observed inner Galactic pulsars may be 
due to the high electron density in this region of the Galaxy 
and therefore larger scattering of the pulse signal, however, 
the latest observations do suggest an intrinsic scarcity of 
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Figure 5. Left panels: theoretical radial distributions of pulsars for a range of models. All model radial distributions are normalised to 
Model A. Right panels: theoretical space velocity distributions for the same range of models. The velocity distributions depicted here 
represent the pulsar population space velocities with respect to the Galactic centre. The solid line in the bottom-right panel depicts the 
3D space velocity distribution of young pulsars derived from observations (Hobbs et al. 2005). All velocity distributions are normalised 
to unity. 



central pulsars (Lorimer et al. 2006). Model A provides a 
good comparison to the observed radial pulsar distribution 
for the inner regions of the Galaxy but has too many pul- 
sars and is too extended beyond ~ 4 kpc. In terms of shape, 
Model C best represents the observations. However, Model 
C peaks further from the Galactic centre by 1 — 1.5 kpc. This 
suggests that the ideal initial distribution would of the form 
derived by Yusifov & Kucuk (2004) from observations of OB 
stars but scaled so that the distribution peaked at a radius 
of ~ 4 kpc. 

The scale heights for Models A, B and C can be com- 
pared in Table[2l We see that Model B has systematically the 



largest scale heights while Model C has the smallest. How- 
ever, the trends observed for Model A - the relative scale 
heights of the various pulsar populations - are consistent 
across the models 

We next demonstrate what effect modifying the Galac- 
tic region of interest has on Model C by considering pulsars 
only out to \z\ = 2 kpc in Model C2 and out to \z\ = 20 kpc 
for Model C3 (as opposed to \z\ = 10 kpc for Model C). We 
find that reducing the height of our 'Galaxy' by a factor of 
five approximately halves the calculated scale heights for all 
pulsar populations. On the other hand, doubling the height 
considered does not significantly affect the calculated scale 
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Table 3. Observed scale heights (in kpc) for different types of 
pulsars. We have not taken account of selection effects when cal- 
culating these values. Standard pulsars are those pulsars which 
satisfy Equation 1341 while MSPs are those pulsars which have 
spin-periods P ^ 0.02 s. We note that there is uncertainty in 
some of these numbers owing to small number statistics and the 
clumpy distribution of the pulsars. For example, there are only 19 
isolated MSPs and the difference in height between the 12th and 
13th most distant (in terms of height from the plane) is 0.07 kpc 
while the average distance between the first twelve is 0.02 kpc. 



Type 


All 


Standard 


MSP 


Both 


0.40 


0.39 


0.41 


Isolated 


0.41 


0.39 


0.26 


Binary 


0.38 


0.17 


0.48 



heights (except perhaps the isolated MSP population which 
suffers from small number statistics). However, it does not 
appear to greatly affect the relative scale heights of pulsar 
sub-populations and certainly does not switch any trends 
noted in the models. Beyond this limit only the results of 
highly energetic systems may still vary. Therefore, factors 
which limit the region of the Galaxy observed (or consid- 
ered), such as the numerous selection effects which occur in 
radio pulsar observations, can modify the underlying pulsar 
population scale heights within \z\ < 10 kpc of the Galac- 
tic plane (as discussed by many works including Taylor & 
Manchester 1977 and Narayan & Ostriker 1990). We note 
that in terms of Galactic pulsar observations there is the 
limit of ~ 1.75 kpc beyond which dispersion measure dis- 
tance estimates of pulsars break down (see Manchester et 
al. 2005). 

We now have Model C, a suitable model for which we 
may compare pulsar scale heights to observations (the lat- 
ter values are given in Table [Sjl . Model C2 is roughly con- 
sistent with the observed scale heights, although on aver- 
age the model values are greater than the observed values. 
The trends when considering all pulsars are similar but this 
breaks down for the MSP population where the model pre- 
dicts a greater scale height for isolated MSPs than for binary 
MSPs which is opposite to the observed MSP scale heights 
(although, see Section I4.4.ip . This is, however, consistent 
with our previous simple analysis in Section [T] from the bi- 
nary disruption formation mechanism of isolated MSPs. An- 
other difference between Model C2 and observations is the 
relative number of isolated to binary MSP systems, in Model 
C2 ~ 99% of MSPs are found within a binary system while 
a direct observational comparison shows ~ 70% of Galactic 
disk MSPs in binary systems. These two differences between 
our model and observations indicate that our mechanism for 
producing isolated MSPs - binary disruption in a SN event 
- can not be the sole (or even dominant) production mecha- 
nism. We explore this line of thought further in Section [4.41 

4.2 Model Galactic gravitational potentials 

We now explore what effect modifying the assumed Galac- 
tic gravitational potential has on the pulsar scale heights 
and final distributions (radial and space velocity). Model 
D assumes a potential of the form described by the KG89 
model (Equation [21} . The scale heights of Model D are all 
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Figure 6. Theoretical radial distributions of pulsars for a select 
few models all normalised to unity. Here we restrict ourselves to 
only consider those pulsars with P > 10~^^ s s~^ which com- 
pares to the observed sample used by Yusifov & Kucuk (2004: 
solid line). Also provided is the initial radial distribution given in 
Paczynski (1990: dotted line), which is used in Models A, D and 
E, and the initial radial distribution of Yusifov Sz Kucuk (2004: 
dash dotted line) used in Model C. 



slightly less than their Model A counterparts which used 
the P90 model. This suggests that over time stellar systems 
may diffuse less efficiently in Model D than in Model A. Fig- 
ure 13] gives some indication of the cause of this difference. 
The lower panel of provides the model Galactic gravitational 
force towards the plane with respect to the height above 
the plane. We see that above a height of |2;| ~ 1 kpc the 
KG89 model has a slightly greater attractive force than the 
P90 model. However, for the inner regions the situation is 
reversed and indeed if we calculate the scale heights for pul- 
sars with 1^1 < 1 kpc, we find that the scale height behaviour 
for Model D relative to Model A is also reversed. As shown 
in Figures [5] and [6] the radial distribution does not differ 
greatly between Models A and D: small differences to note 
are that the distribution of Model D decays more rapidly 
than Model A while Model D contains less pulsar systems 
within 10 kpc of the Galactic plane. The velocity curve of 
Model D does not change greatly from Model A as there is 
only a small increase in systems with space velocities at the 
lower end of the distribution, perhaps reflecting the rapid 
decay in circular velocity of Model D with respect to Model 
A (see Figure |4}. 

Model E assumes the Galactic gravitational potential 
of Xue08 (Equations 1281 to l30p . This model shows a large in- 
crease in scale height values compared to both Models A and 
D. The reason for this can once again be seen from Figure [4] 
which shows that the Xue08 Galactic gravitational model 
can not retard the movement of systems out of the Galactic 
plane as effectively as the P90 (Model A) or KG89 (Model 
D) models. Figure [5] shows that the number of systems re- 
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tained by Model E within \z\ ^ 10 kpc is less than within 
both Models A and D. The 3D space velocity curves of Mod- 
els A, D and E all peak at roughly the same value but the 
Model E distribution possesses a much steeper slope either 
side of its peak value. The flatness of the rotation curve 
assumed in Model E, depicted in Figure ID causes such a 
narrow distribution in velocity space. Model E highlights 
the importance of the assumed Galactic gravitational po- 
tential in modelling Galactic population kinematics (as also 
discussed by Kuijken & Gilmore 1989, Dehnen & Binney 
1998 and Sun & Han 2004). 

4.3 Kicks 

Considering the uncertainty involved in the true form of the 
distribution of kick speeds given to NSs at birth (as men- 
tioned in Section 2.2.1) we next investigate how changing 
the dispersion of the assumed Maxwellian distribution af- 
fects our results. Recalling that Model A used Va = 190 km 
s~^ we first compare with Model F which uses Va = 550 km 
s~^ as an extreme illustration. We see from Table[2]that the 
scale heights in Model F are much greater than in Model 
A, with increases by as much as a factor of 2. While it is 
expected that pulsars can move further from the galactic 
plane in Model F it also means that more objects escape 
from the Galaxy: the number of pulsars in Model F within 
10 kpc of the plane decreases by a factor of 10 compared to 
Model A. This decrease in numbers can skew the expected 
outcomes and is evident when looking at the radial distri- 
butions in Figure \S\ For example, less binary systems are 
kicked out of the \z\ ^ 10 kpc Galactic region than isolated 
pulsars - owing to binaries being heavier on average and the 
binary orbit absorbing a fraction of the energy injected by 
the kick - so we find in Model F that binary pulsars have 
a greater scale height than isolated pulsars (the opposite to 
Model A). This is not the case for MSPs although the dif- 
ference between the isolated and binary MSP scale heights 
has decreased compared to Model A. We note that the ratio 
of standard isolated pulsars to standard binary pulsars is an 
order of magnitude greater for Model F than for Model A 
owing to a greater incidence of binary disruptions in Model 
F. Therefore, as established in previous works (e.g. Taylor & 
Manchester 1977; Hills 1983) , the assumed SN kick velocity 
is an extremely important factor, especially when comparing 
model pulsar kinematics to observations. 

When examining the resultant radial distribution in 
Figure [5] we find that Model F is much broader than Model 
A - an intuitive results owing to the increased distance that 
the pulsars move - and certainly Model F is not a good 
representation of the observed distribution (not that it was 
expected to be). Comparing the space velocity distributions 
of Models A and F in Figure [5] it is surprising to see the 
relatively high number of pulsars in Model F that are trav- 
eling at the distribution peak speed (~ 250 km s~^). This 
population with similar velocities includes isolated and bi- 
nary pulsars. We remind the reader that when discussing 
the velocity distribution here we mean the actual velocity 
each system has with respect to the Galactic centre - we do 
not account for the local standard of rest (solar motion). 

In Figure [S] we also compare the 3D space velocity dis- 
tributions of Models A and F with the 3D space velocity 
distribution derived by Hobbs et al. (2005) from pulsar ob- 



servations. An important distinction to make is that the 
Hobbs et al. (2005) sample was restricted to pulsars with 
characteristic ages < 3 Myr. Thus, it is intended to be a dis- 
tribution of pulsar birth velocities and was used by Hobbs 
et al. (2005) to suggest that = 265 km s~^ in the SN 
velocity kick distribution. By comparing this with our mod- 
els we can gauge the effect that the Galactic potential and 
binarity have on the form of the pulsar velocity distribution 
as the population evolves. 

Comparing the Model A pulsar velocity distribution (at 
a population age of 10 Gyr) with the Hobbs et al. (2005) 
birth velocity distribution shows changes in the peak (shifted 
to lower velocity in the model) and shape. The shift of the 
peak can be mostly attributed to the difference in the aver- 
age age of the two pulsar populations and the low dispersion 
value assumed in the Model A velocity kick distribution. The 
age difference allows fast moving systems in Model A to have 
time to leave the Galaxy and thus be culled from the final 
velocity distribution. Also, over time, the pulsar velocities 
are retarded by the Galactic potential which shifts the final 
pulsar velocity distribution to lower velocities. In terms of 
shape we find that Model A is a more focussed distribution - 
the model peak is more acute and the distribution as a whole 
is narrower. This difference is likely due, in part, to the bind- 
ing energy of the host binaries impinging on the SN velocity 
kick (as discussed in Hills 1983 and Bailes 1989) and causing 
a greater number of systems to have similar final space ve- 
locities than otherwise. Disruptions triggered primarily by 
mass loss will act to increase the proportion of low veloc- 
ity pulsars while the absorption of large kick velocities by 
the binary binding energy may also skew the distribution to 
smaller final pulsar velocities. This narrowing of the model 
pulsar velocity distribution compared to the observed distri- 
bution has been found in other works, most notably that of 
Dewey & Cordes (1987). They attributed the difference to 
errors in pulsar distance measurements, which will broaden 
the distribution, and that non-Maxwellian processes may 
be more important in producing pulsar velocities than their 
models assume (that is, nascent NS receiving a kick selected 
from a Maxwellian distribution with < V >= 90 km s~^). 
It appears that the difference between the velocity distribu- 
tion of models and observations results from a combination 
of the selected pulsar population used to derive the observed 
pulsar velocity distribution (see Hobbs et al. 2005), errors in 
pulsar velocity and distance measurements, and the binding 
energy of host binary systems affecting the resultant pulsar 
run-away velocity. 

To remove the binary orbit effect and highlight the ef- 
fect of age evolution on the pulsar velocity distribution we 
have created Model G which evolves a population of single 
stars according to the Galactic setup described for Model 
A but with V<T = 265 km s~\ With every system evolved 
within Model G being isolated from birth we are now able 
to compare the resultant velocity distribution of a popula- 
tion of pulsars which receive uninhibited SN kick velocities 
drawn directly from the suggested Hobbs et al. (2005) SN 
kick distribution. We now see in Figure [5] that the distribu- 
tion closely resembles the Hobbs et al. (2005) distribution 
in shape but is shifted to lower velocities. The final pulsar 
distribution is best fit by a Maxwellian distribution with 
Va = 140 km s-^ 
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Table 4. Scale heights (in kpc) of MSPs and their companions 
for Model C (see Section 4.4.1) and its variants (see Sections 4.4.4 
and 4.4.5). 



Model 


C 


c' 


c" 


MSP-MS 


0.74 


0.75 


0.74 


MSP-WD 


0.83 


0.84 


0.86 


MSP-NS 


1.96 


1.93 


2.06 


MSP-BH 


0.11 


0.13 


0.10 


Isolated MSP 


1.76 


1.76 


1.61 


Ablation 


C 


c' 


c" 


Binary MSP 
Isolated MSP 


0.87 
0.75 


0.87 
0.75 


0.91 
0.77 



4.4 Millisecond pulsars 

In Paper I we were primarily interested in the production 
of MSPs within the P-P diagram. We now continue our 
exploration of the MSP population by examining in more 
detail the Galactic MSP distributions and in particular fo- 
cussing on the behaviour of isolated MSPs and those with 
MS star, WD, NS or BH companions. In doing this we fo- 
cus solely on Model C. To begin, we extend our evaluation 
of scale heights in Table [2] with those of the MSP popula- 
tions (given in Tableland discussed in Section r4.4.ip . The 
scale heights in Table |4] are supplemented by Figure [7] which 
provides the scale height for each MSP population as a func- 
tion of Galactocentric radius. Also shown is the Galactic x 
and z parameter space of MSPs: for all MSPs (see Figure [8]) 
and those that reside above a magnetic field cut-off (see Fig- 
ureO . The population of MSP-BH binaries is then discussed 
in detail within Section l4A2] We then look at the MSP pop- 
ulation recoil velocities and space velocities in Section [4.4.31 
and make some limited comparisons to previous work and 
observations. To further our investigation into how different 
model assumptions affect our pulsar population kinematics, 
we modify Model C, our favoured model thus far, to account 
for: a greater birth |2maa:i| range (Model C in Section r4.4.4|) : 
a greater age of the Galaxy (Model C in Section r4.4.5p : and 
a higher resolution sample (in Section [4. 4. 6p . 

4-i-l Model C MSP scale heights and Galactic spatial 
properties 

Looking at the scale height values in Table |4] for Model C we 
see that as expected the binary MSP population with the 
greatest scale height is the MSP-NS systems, in which two 
SNe kicks occur. These double compact systems, however, 
are much rarer than the MSP-MS or MSP-WD systems and 
therefore the results are less statistically significant. The rel- 
ative numbers of MSP-NSs compared to MSP-WD systems 
(the largest MSP population) is 0.003. For MSP-MS systems 
(the second most numerous MSP population) the relative 
number is 0.044 per MSP-WD system. We find similar scale 
heights for the MSP-MS and MSP-WD systems although 
the former are systematically smaller owing to the popu- 
lation being younger on average. Recently there have been 
suggestions that asymmetric mass-loss during the asymp- 
totic giant branch phase (Spruit 1998) give rise to WD re- 
coil velocities of the order of a few kms~^ (Fellhauer et al. 



2003). Such kick velocities have been raised a possible ex- 
planations of the apparent deficit of WDs in open clusters 
(Fellhauer et al. 2003) and the radial distributions of WDs in 
globular clusters (Heyl 2007; Davis et al. 2008). Currently 
we do not include this possibility in our models but note 
that it would presumably lead to a modest increase in the 
MSP-WD population scale height. 

MSP-BH systems are found to have a small Galactic 
scale height. This arises due to the orbital parameters re- 
quired in order to form these systems which we examine in 
further detail below (see Section |4]4]2}. Also, we remind the 
reader that we currently assume BHs do not receive kicks 
during their formation. As shown in Tables [5] and [3] the iso- 
lated MSP population has a scale height of 1.76 kpc. These 
MSPs emerge from disrupted binary systems and although 
the kick at the time of disruption may be large it is not the 
MSP which is exploding at that point. Therefore the MSP 
is considered by our kick routine to be the secondary star 
which, as shown in Section [2.2.31 receives (on average) only 
a small increase in momentum. This results in the lower 
scale height of isolated MSPs compared to MSP-NS binary 
systems (albeit only slightly less than the MSP-NS value). 
Furthermore, we note that for the binary system to survive 
the first SNe, allowing mass transfer onto the progenitor 
MSP, the resultant velocity kick at this point must be rel- 
atively small (we find Vkick of approximately 80 km s~^ or 
less). This is in accordance with many other population syn- 
thesis works, including StoUman & van den Heuvel (1986), 
Iben, Tutukov & Yungleson (1995) and Ramachandran & 
Bhattacharya (1997). 

Previous results shown in Section 14.11 placed doubt on 
isolated MSPs formed via the disruption of binary systems 
being the sole 'type' of isolated MSPs - there must be an- 
other formation mechanism. One such mechanism that ex- 
ists in the literature is the ablation model (Eichler & Amir 
1988; Ruderman et al. 1989) based on observations such as 
those of van Paradijs et al. (1988). Here the assumption is 
that the MSP is produced as a result of mass-transfer from 
a MS companion in what would be a low-mass X-ray bi- 
nary. Then at some point the mass of the MS star becomes 
low enough that it is destroyed, or ablated, by the highly 
energetic radiation flowing from the rapidly spinning pulsar 
(van Paradijs et al. 1988; Tavani 1992). We calculate that the 
timescale for the destruction of the MS companion star in 
this manner should take of order ~ 5 Myrs once the compan- 
ion is below a mass of ~ 0.02 M0 (see Appendix^). Thus 
we propose a simple model to belatedly estimate the impact 
of ablation on our results where we assume that any MSP 
with a MS companion of mass less than O.OOSMq (to be on 
the safe side) is in fact an isolated MSP. With the inclusion 
of ablation we find that the percentage of isolated MSPs 
increases from 1% to 36%. This new value is in rough agree- 
ment with observations where it is estimated that one third 
of the MSPs are isolated (Huang & Becker 2007). Iben, Tu- 
tukov & Tungleson (1995) similarly found good agreement 
with observations for binary to isolated ratios when assum- 
ing ablation of MSP companions. We see from the last two 
rows in Table |4] that the isolated and binary MSP popula- 
tions now have comparable scale height values (in fact the 
isolated scale height is now slightly the lower of the two). 
Therefore the kinematics of the binary and isolated MSP 
populations are now similar. This last point is actually con- 
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Figure 7. The radial variation of pulsar scale heights for a range 
of MSP populations in Model C. 



sistent with the observations of MSPs, which via statistical 
arguments show no difference in binary and isolated MSP 
kinematics (Lorimer et al. 2007). From this simple test we 
see that the low mass companions to MSPs do occur and 
that the ablation process deserves serious consideration in 
future models. 

In Figure [7] we show how the scale heights of the MSP 
populations vary with Galactocentric radius. We note that 
the region of the Galaxy where the populations are most nu- 
merous is between 4 — 6 kpc from the Galactic centre. The 
top panel of Figure [7] depicts the similarity of MSP-MS and 
MSP-WD kinematics. It is only out beyond ~ 13 kpc that 
the two populations diverge, and this is only due to low 
number statistics which begin to plague the MSP-MS re- 
sults. Low number statistics have a much greater influence 
in the middle panel of Figure [T] For example, the highest 
number of systems in a radial bin (1 kpc in width) for the 
MSP-NS population contains 39 systems while the lowest 
only 3. The MSP-NS and isolated MSP systems have similar 
scale heights throughout the majority of the Galaxy (after 
accounting for statistical uncertainty) and systems can be 
found far from the plane. The MSP-BH systems on the other 
hand are all found close to the Galactic plane. When ac- 
counting for the ablation of MSP companions we flnd a very 
similar distribution of binary and isolated MSPs throughout 
the entire Galaxy. 

It is also interesting to compare the spatial Galactic 
X — z distributions of the MSP populations. This is shown 
in Figure [8] for the Galactic a;2-plane and emphasises what 
we have already seen in Tableland Figure [7] isolated MSPs 
and MSP-NS binaries have quite extended distributions (rel- 
ative to their numbers) , MSP-BH systems reside close to the 
plane and the majority of MSPs are found with WD com- 
panions (MSP population numbers relative to MSP-WD sys- 
tems are given in the lower right corner of each panel) . What 
is surprising in Figure|8]is the large number of MSP-WD sys- 
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Figure 8. Galactic x and z coordinates for Model C MSP systems 
within a radius R = i^/a;^ + j/^ < 30 kpc of the Galactic centre 
and \z\ ^10 kpc from the Galactic plane. The bottom left corner 
of each panel gives the MSP companion type and relative number 
of that system compared to the MSP-WD systems. Due to their 
lack of numbers the points for MSP-NSs, MSP-BHs and isolated 
MSPs are larger than for the MSP-MS and MSP-WD systems. 



terns out to jzj = 10 kpc. This suggests that there may be 
many MSP-WD systems lost from - but surrounding - the 
Galaxy. We next investigate the result of imposing a limited 
selection effect on the MSP population where we only con- 
sider pulsars that have B > 6 x 10^ G. This magnetic field 
value is a suggested limit (Zhang & Kojima 2006) of the 
required field strength to turn on (or off) the pulse mecha- 
nism (see Paper I). Figure |9] shows the field strength limited 
MSP population and the result in comparison to Figure |8] is 
dramatic. The entire MSP-BH population now disappears, 
which is also almost the case for the MSP-MS population 
where only two systems are left. The relative numbers of 
both isolated MSPs and MSP-NSs have now increased com- 
pared to the MSP-WD systems (see values on Figures [8] and 
[9|. Clearly many MSPs in Model C accrete enough mass 
to cause a large decay in the magnetic field. In particular 
every pulsar within a MSP-BH system has accreted more 
than ~ 0.04 AIq of material which is the typical amount of 
mass it takes for our model pulsar magnetic fields to decay 
below B — Q X 10^ G (see Paper I). This is compared to 
other works which assume AM > 0.1 Mq is required for 
MSP production (e.g. Willems & Kolb 2005). 



4.4.2 Model G MSPs and BHs 

Although other works have detailed BH and pulsar binary 
evolution in varying detail (Narayan, Piran & Shemi 1991; 
Lipunov et al. 1994; Pfahl, Podsiadlowski & Rappaport 
2005; Lipunov, Bogomazov & Abubekerov 2005) we eval- 
uate the accretion history of MSP-BHs and why these sys- 
tems reside close to the Galactic plane. To place this into 
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Figure 9. As for Figure |8] but now restricted to include only 
pulsars that satisfy B 6 X 10*^ G. 

context we explore the initial orbital period and initial pri- 
mary mass (the MSP progenitor and initially the more mas- 
sive star) parameter space in Figure 1101 This figure is also 
designed to show which systems reside in a range of ob- 
servable orbital periods (that is orbital periods which would 
be observed now, if the age of the Galaxy is 10 Gyr). Fig- 
ure [TO] also gives the secondary mass (BH progenitor and 
initially the less massive star) range for each binary sys- 
tem depicted. What we find in Figure [10] is three distinct 
groups in both initial orbital period and final orbital period. 
Each of the three regions have different evolutionary path- 
ways which leave the systems close to the Galactic plane. 
We now describe these evolutionary pathways. Firstly, how- 
ever, we note the healthy number of MSP-BH systems (few 
hundred) produced within our model. This is in contrast to 
Pfahl, Podsiadlowski & Rappaport (2005) who suggest that 
BHs with recycled pulsar companions are rare. The primary 
differences between our model and that of Pfahl et al. (2005) 
- in terms of BH formation - is their assumption that SN 
kicks are given to BHs, the assumed evolutionary parameter 
values within the CE phase and massive star wind mass-loss 
and accretion. In varying the CE parameters Pfahl et al. 
(2005) show that ~ 100 MSP-BH systems can be produced 
(with extreme CE parameter choices), however, a favoured 
model estimates only 5—10 MSP-BH systems within the 
Galaxy. This compares well with the number of systems we 
find with orbital periods less than ~ 10 hr, which is roughly 
the maximum orbital period limit Pfahl et al. (2005) find for 
MSP-BH systems. We also include an evolutionary pathway 
which forms an MSP-BH system that does not include a CE 
phase, something not considered by Pfahl et al. (2005; who 
do, however, comment on this scenario). 

Those systems which begin their lives with initial or- 
bital periods less than 10 days all have an initial primary 
mass of Mii > 40 Mq and an initial secondary mass of 
M21 > 20 Mq. These systems end with the largest BH 



masses (around ~ 13 Mq) of all the MSP-BH populations 
and their orbital periods are grouped near 1000 days. These 
systems were not found by Pfahl et al. (2005) possibly ow- 
ing to the inclusion of SN kicks during BH formation. The 
general evolution pathway of these systems goes as follows. 
The initial orbital separation is of the order of 80 -Rq or less 
and the massive primary star evolves to fill its Roche-lobe 
within a few Myr. This leads to a phase of steady mass 
transfer lasting 1 — 2 Myr and ending with the primary as 
a naked helium star with a mass of about 10 Mq . During 
the phase the secondary accretes approximately 80% of the 
transferred material with the remainder lost from the sys- 
tem. The orbital separation at this point is typically 200 Rq 
and subsequently increases further owing to winds from the 
helium star and the now massive secondary. At a system 
time of ~ 5 Myr the primary undergoes a SN explosion and 
becomes a NS. We find that velocity kick magnitudes of 
Vkick ^ 80 km s~^ allow the system to remain bound. Be- 
yond the first SN the secondary evolves quickly and loses a 
large proportion of its matter in a wind, some of which is 
accreted by the NS companion. The secondary evolves via a 
naked helium star phase to explode as a SN and leave a BH 
remnant. At this point we have an eccentric MSP-BH system 
which has received one mild SN velocity kick in its lifetime 
and has typical component masses of 2 and 10 — 13 Mq, for 
the NS and BH respectively. The orbital separation is in the 
range of 1000 — 4000 -Rq (depending on the precise details 
of the kick velocity and the mass-loss history) . 

Those MSP-BH systems with 10 < Porbi < 100 days as 
seen in Figure \TU\ end their lives with a large range of BH 
masses extending from 3 Mq through to 11 Mq in tight 
orbits around their MSP companion (Porbf < 20 days). 
It is this population of MSP-BHs which are most likely 
to coalesce at and around the age of the Galaxy (simi- 
lar to Pfahl et al. 2005). Initial primary masses are in the 
18 < Mu/Mq < 30 range and secondary masses are typ- 
ically 10 < M2i/A'/Q < 20. The initial orbital separation 
ranges from 100 — 300-Rq. Early evolution proceeds simi- 
larly to that of the previous group: non-conservative mass 
transfer from the primary to the secondary accompanied by 
an increase in the orbital separation and ending with the pri- 
mary as a naked helium star. The primary then undergoes 
a SN and becomes a NS at a system time of about 8 Myr. 
We find that generally these systems can survive slightly 
larger SN velocity kicks than the systems described in the 
previous group. The companion is now a massive MS star 
(~ 30 Mq) and subsequently fills its Roche- lobe while cross- 
ing the Hertzsprung Gap. This initiates dynamical-timescale 
mass transfer leading to a common-envelope phase and the 
creation of a tight binary comprised of the NS primary 
(~ 2AfQ) and a naked helium star secondary (~ 10 Mq). 
We note that systems in the first group avoid this second 
Roche-lobe filling event because the secondary is more mas- 
sive and loses mass in a wind at a greater rate leading to 
more substantial orbit expansion after NS formation. After 
emerging from the common-envelope the NS then accretes 
material from the wind of the companion to become a MSP. 
This ends when the companion becomes a BH. The final 
MSP-BH binary will have an orbital separation of less than 
10 -Rq and systems such as this may coalesce within a Hub- 
ble time. 

The MSP-BH systems with Porbi > 100 days end with 
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Figure 10. MSP-BH initial population parameter space of Model 
C, in particular the initial orbital period and zero-age main- 
sequence primary star (MSP progenitor) mass. Provided are 
ranges of the initial companion mass (BH progenitor). Those sys- 
tems with initial orbital periods, Po^bi > 100 days end with final 
orbital periods P^-ch > 10000 days. Those with Porbi < 10 days 
end with 20 < Porb < 10000 days and those with 10 < Porbi < 
100 days end with Po^b < 20 days. 
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Figure 11. Model C MSP-BH pulsar \z \ — Porbf parameter space. 
There are three distinct groups in final orbital period. The group 
on the left corresponds to systems with 10 < Porbi < 100 days, 
the central group corresponds to systems with Porbi < 10 days 
and the right hand group corresponds to systems with Porbi > 
100 days (see Figure [To)l . 



orbital periods of 1000 days or more (see Figure [TT] for the 
final orbital period range). We note that the smallest pri- 
mary and secondary masses belong to this group. Once again 
mass-transfer occurs prior to the first SN event but owing 
to the wider orbit this is initiated later (~ 15 Myr) than in 
the previous cases and when the primary is a giant star. The 
orbital separation when the primary undergoes a SN (to be- 
come a NS) is typically 2000 — 3000 Rq which means that 
relatively smaller kicks are required if the system is to re- 
main bound and proceed to become a MSP-BH binary. We 
find that kicks of the order of 20 km s~^ or less are necessary 
(slightly larger if the kick is well directed). After NS forma- 
tion the secondary is a MS star with mass of approximately 
20 Mq . The secondary then evolves off the MS and transfers 
some material to the NS before ending its life as a BH of 
mass less than 5 Mq . 

The above analysis shows that the most likely MSP- 
BH systems to be created are those in which the first SN 
- the only one assumed to impart a velocity kick onto the 
compact remnant - produces a small velocity kick, which 
is why these systems are found to hug the Galactic plane 
as suggested by Narayan et al. (1991). In fact, compared to 
the other MSP binary populations the MSP-BH systems ef- 
fectively represent a different kick distribution, in that the 
distribution of kicks given to systems that remain bound is 
distinct. As touched on in the evolutionary descriptions this 
is also true internal to the MSP-BH population, where the 
effective kick distribution for systems that remain bound is 
different for each of the three period groupings we identi- 
fied in Figure 1101 This is depicted indirectly in Figure 1111 
Here we see the MSP-BH height from the Galactic plane 
and the populations are designated by their grouping in the 
final orbital period parameter space. Each population has a 
different scatter in \z\, which can be used as an indicator for 
the average strength of the SN velocity kick. We see that 
the majority of those small orbital period MSP-BH systems 
are further off the plane than the extremely long period 
MSP-BH systems, suggesting that as expected from Bailes 
(1989), the close binary systems can survive larger kick ve- 
locities than the larger binary systems (which was outlined 
in the evolutionary examples). Only three MSP-BH systems 
are found beyond 10 kpc from the Galactic plane. 

The MSP-BH orbital period distributions as shown in 
Figure [TT] are remarkably distinct and perhaps surprisingly 
not smeared out by our use of random birth ages. This is due 
to the vast orbital period differences between these popula- 
tions and the time scales these populations evolve on. The 
orbit of MSP-BH binary systems, after the formation of the 
BH, can only shrink in time owing to gravitational radia- 
tion (Landau & Lifshitz 1951; Hulse & Taylor 1985; Hurley, 
Tout & Pols 2002). However, the time-scale on which this 
decrease occurs is greatly dependent on the size and eccen- 
tricity of the orbit. Long period binary systems have very 
large timescales for orbital parameter change and thus re- 
main as long period systems over a Hubble time. The very 
close systems (separation ^ 10 Rq) will shrink more rapidly 
and may even coalesce within a Hubble time. Therefore, the 
long period systems stay long and the short period systems 
only get shorter and as a result the MSP-BH systems stay 
within their orbital period groups as they evolve throughout 
the Galaxy. Thus we observe three distinct MSP-BH groups, 
a result differing some what from the orbital period distri- 
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bution of Pfahl et al. (2005) who find that most MSP-BHs 
have orbital periods of 1 — 6 hr. 

4.4-3 Model C MSP recoil and 3D space velocities 

We now examine the MSP population in velocity and or- 
bital period parameter space. In terms of velocity we con- 
sider both the recoil velocity and the space velocity of the 
binary systems. The recoil velocity is defined as the change 
in velocity of the binary centre-of-mass owing to the SN 
explosion that created the NS (that went on to become the 
MSP). The space velocity is the velocity of the binary within 
the Galaxy at the time when the Galactic age is 10 Gyr. In 
calculating the final space velocities the solar motion around 
the Galactic centre is accounted for by removing the local 
standard of rest (LSR) velocity of ~ 220 km s~^ (Dehnen 
& Binney 1998). The orbital period is taken as the final or- 
bital period at a Galactic age of 10 Gyr. In a similar vein 
to Section [4AT] we examine the parameter space when con- 
sidering all pulsars and then examine it again after limiting 
the sample population to MSPs that have B ^ 6 x 10^ G. 
The results are shown in Figure [121 

The recoil velocities for all MSP binaries are shown in 
Figure [T2k ) . The first item to note is the MSP-BH systems 
which all have low recoil velocities but cover a large range 
of final orbital periods. Such a distribution is not surprising 
given our detailed analysis of such systems in Section [4.4.21 
Also not surprising is the rather large recoil velocity range 
of double-NS systems. The typical total recoil velocity inci- 
dent on such systems is greater than 200 km s~^. We can 
also see from Figure I12h 'l that these systems are likely to 
be eccentric rather than circular. For MSP-NS systems that 
receive large recoil velocities (> 450 km s~^) there appears 
to be a lower limit to the possible final orbital period. The 
initial orbital period of these systems is very important in 
determining the evolution outcomes and the appearance of 
the final parameter space (Tauris & Bailes 1996). Also, in 
a related manner and as discussed for MSP-BH systems in 
Section 14.4.21 the onset of mass transfer and the details of 
the common-envelope phase are crucial factors. What we 
find is that a significant proportion of the double-NS popu- 
lation end up with extremely small periods (and a range of 
eccentricities) and coalesce rapidly (within a few Myr after 
double-NS formation) similar to that found by Belczynski, 
Bulik & Rudak (2002). This leads to the orbital period gap 
observed in Figure fT^ ). We leave further discussion on these 
systems for future work (Kiel, Hurley & Bailes, submitted). 
Turning to the MSP- WD systems we see that these typically 
receive rather low recoil velocities with the average value 
being less than 100 km s~^ (much less than \/^Vcr) and in 
accordance to previous population synthesis results of Ra- 
machandran & Bhattacharya (1997), Phinney & Kulkarni 
(1994), Lyne et al. (1998) and Sun & Han (2004). We also 
see that a similar but opposite trend occurs for MSP-WDs 
as did for the MSP-NSs in that for large recoil velocity val- 
ues there appears to be an upper limit to the possible final 
orbital period. Again this is related to the orbital evolution 
and in particular whether a system enters common-envelope 
evolution (and survives without coalescence) or not. 

Figure 112b ) shows the recoil velocity and final orbital 
period parameter space for the magnetic-field limited MSP 
population. We see that the population has been signifi- 



cantly thinned out. In particular the entire MSP-BH popu- 
lation has been removed as have the low-period MSP- WD 
systems. It is also possible to compare our findings to the 
results of Tauris & Bailes (1996: see their Figure 2c) who 
followed the formation of MSPs using stellar and binary evo- 
lution algorithms that were quite advanced for their time. 
Compared to Tauris & Bailes (1996) we find less systems 
with orbital periods greater than a day. However, for the 
MSP- WD population we observe a similar trend of orbital 
period to recoil velocity: the smaller the orbital period the 
greater the range in possible recoil velocity of the system. 



In Figure [T2b l we look at the final space velocities and 
orbital periods for all MSP binaries. There is much similar- 
ity between the velocities given to each system (their recoil 
velocities) and their LSR Galactic motion. The form of this 
parameter space is therefore governed by the same evolu- 
tionary phases that dictated the appearance of Figure [12^). 
In Figure I12t l) we show the space velocity-orbital period 
parameter space distribution of MSPs with B ^ 6 x 10^ G. 
Included for comparison are pulsar proper motion observa- 
tions which, convolved with distance estimates, give rise to 
observed transverse velocities. From the ATNF pulsar cata- 
logue (Manchester, Hobbs, Teoh & Hobbs 2005) there are at 
present 28 Galactic disk pulsars with spin periods less than 
0.02 s that have measured orbital periods and estimated 
transverse velocities. Although we do not directly compare 
total model space velocities to observed pulsar proper mo- 
tions, useful information can still be gleaned from simple 
comparisons between the two noting that the transverse ve- 
locities are a lower limit to the true space motion (although 
measurement errors not included within Figure 1121 espe- 
cially in distance calculations, cloud this picture slightly). 
Firstly, once accounting for the LSR, we can see that many 
of the model MSP systems travel with speeds within the 
typical stellar velocity range of approximately ±16 km s~^ 
(as given by Dehnen & Binney 1998). However, there ap- 
pears to be an overabundance of model MSP binaries at low 
velocities. The model also fails to produce enough of the 
fast moving MSPs with large orbital periods. One reason 
for this may be the Vcr = 190 km s~^ assumed in Model 
C which is lower than the value suggested from observa- 
tions {Va = 265 km s"^: Hobbs et al. 2005). On the other 
hand, Figure [T2Ji) shows a large range of space velocities for 
MSP-NS systems. Perhaps surprisingly some of these sys- 
tems even have velocities close to the LSR, although not so 
surprising according to Dewi, Podsiadlowski & Pols (2005) 
who suggest that DNSs receive small kicks. This is of par- 
ticular interest for understanding the nature of the double 
pulsar J0737-3039 (Burgay et al. 2003) which is observed 
to have a transverse velocity of 30 km s^^ or less (Kramer 
et al. 2006). Considering that the system will have experi- 
enced two supernova events this has been taken as evidence 
for little or no velocity kicks within this system. However, 
Kalogera et al. (2007) have described models which show 
that kick velocities of 100 km s~^ or more are still possi- 
ble. Our results agree with this in that it is not necessary 
to make any unusual assumptions regarding kicks in binary 
systems to explain the observed velocities of systems such 
as J0737-3039 (Deller, Bailes & Tingay 2009). 
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Figure 12. Parameter space of velocities (recoil and space) and orbital period for MSPs in Model C. The upper panels show all systems 
while the lower panels are restricted to include only pulsars with _B ^ 6 X 10^ G. The recoil velocities are with respect to the pre-NS 
binary CofM. The space velocities are with respect to the local standard of rest velocity (of ~ 220 km s~^: Dehnen & Binney 1998). 
Identified are MSPs with MS star (0 symbols), giant star (solid squares), WD (plusses), NS (triangles) and BH (stars) companions. The 
large darker points represent those systems with eccentricities greater than 0.1. The filled circles within panel d) represent the 28 MSPs 
with observed transverse velocities and orbital periods taken from the ATNF Pulsar Catalogue. 



4-4-4 Effects of the assumed initial scale height 

Up until now we have set a maximum of l^maxil = 75 pc 
to the initial birth height distribution of binaries, effectively 
modelling a thin disk. We now examine the effect this has 
on the scale height calculations by extending it to |zmaxi| = 
150 pc in Model C . The results are compared to Model 
C in Table ID We find that there is no significant change 
in the calculated scale heights. This agrees with previous 
works, such as Paczynski (1990) or Sun & Han (2004), who 
have suggested that for such kinematically active systems as 



pulsars the initial height above the plane does not greatly 
affect the outcome. 

4-4-5 Effects of the assumed Galactic age 

The age of the Galaxy is an important assumption, espe- 
cially when populations of systems with large differences in 
life times are modelled together. To address this we have 
Model C which assumes the Galactic age is 15 Gyr rather 
than 10 Gyr for Model C. It is the isolated MSPs and MSP- 
NSs whose scale heights change the most appreciably in this 
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new model compared to Model C. This is primarily owing 
to low number statistics (see Figure |8] for the relative num- 
bers of the populations). Otherwise it does not appear that 
the increase in Galaxy age has a significant efi'ect on the 
kinematics of the MSPs. 



4-4-6 Sufficient model resolution? 

Finally, to test whether our previous models have a fine 
enough resolution in the initial parameter space to faithfully 
represent the entire Galactic pulsar population we have ex- 
tended Model C to include 10^ binary systems (a factor of 
100 increase). The only systems for which the scale height 
changed noticeably was the MSP-NS systems, which are rel- 
atively rare and kinematically energetic systems. In all other 
respects it appears that the results for 10^ systems scale reli- 
ably to larger populations. We note here that modelling 10® 
binary systems is equivalent to modelling ~ 10% the mass of 
the Galaxy, assuming binary systems are of interest. To run 
a model of this size takes roughly 4500 CPU hours and even 
when farming the model out to 100 processors (on the Swin- 
burne supercomputei0) it takes almost 2 days to complete. 
Thus it is obviously an advantage when examining a variety 
of evolutionary assumptions one model at a time to be able 
to represent the Galaxy faithfully with fewer systems. 



5 DISCUSSION 

Using our newly developed binkin module for integrating 
the positions of stars and binaries within the Galaxy we 
have worked through a series of models in order to under- 
stand how various options available in the module affect the 
outcomes. This has allowed us to develop a favoured model 
- our Model C. In doing this we have used pulsar popula- 
tions as our yardstick, computing scale heights, radial and 
velocity distributions, and orbital characteristics (in the case 
of binary systems) with limited comparison to observations. 
What we have done is to make predictions in all these areas 
about the particulars of the Galactic pulsar population, as- 
suming that all pulsar systems can be observed. Of course 
this is not the case in reality and our model results cannot 
truly be confronted by observations until we include selec- 
tion effects in our modelling. This will be completed when 
we add our next and final module binsfx to our population 
synthesis code. As such we leave a discussion of the necessary 
selection effects that need to be considered and their treat- 
ment to an upcoming paper focussed on the binsfx module 
(Kiel, Bailes & Hurley, in preparation). This paper will in- 
clude features such as the predicted pulsar P-P diagram for 
distinct regions in the Galaxy (following on from our investi- 
gation of this diagram in Paper I in terms of binary evolution 
parameters). Below we discuss future additions to the bin- 
POP and BINKIN modules in relation to pulsar evolution as 
well as caveats to our current findings. 
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5.1 Accretion induced collapse formation of 
neutron stars 

Further analysis of our MSP populations shows that some 
of the NSs which go on to become MSPs in our models are 
formed from the accretion induced collapse (AIC) of WDs 
(Canal & Schatzman 1976; Nomoto & Kondo 1991). In the 
scenario of Nomoto & Kondo (1991) an O-Ne-Mg WD ac- 
cretes enough material to reach the Chandrasekhar mass, 
the maximum mass possible for a WD to support itself, and 
collapses to form a NS. To date we have allowed NSs formed 
in this way to receive velocity kicks in the same manner as for 
NSs formed in core-collapse SNe. Generally, if the binary sys- 
tem remains bound an AIC NS will continue to accrete ma- 
terial after the SN, a SN which induced an eccentricity into 
the orbit. This formation pathway produces a substantial 
number of MSP- WD and MSP-MS systems with eccentric- 
ities greater than 0.1, that reside within the Galaxy. These 
systems highlight the importance of a correct mass-transfer 
treatment for eccentric binaries (see Paper I and Bonacic- 
Marinovic, Glebbeek & Pols 2008). This is something which 
is not currently accounted for in our models (we make use 
of equations which assume the orbit is circular: see HTP02) 
and it most likely will affect the production and visibility 
of these systems. Of course, if the AIC systems were not 
given any velocity kick, (as has been modelled previously: 
HTP02), or a much lighter kick (as latest models suggest: 
Dessart et al. 2006), then not only would there be many 
more AIC MSP systems but they would all have a greater 
possibility of residing in our Galactic target area and the 
population scale height would be lowered. They would also 
typically have smaller eccentricities. 

Although we do not deal directly with low-mass X-ray 
binaries within this work it is possible for us to compare 
the observed scale heights of such systems - which suffer 
from less selection effects than pulsar observations - with our 
model MSP-MS scale height calculations. This is assuming 
that low-mass X-ray binaries are the progenitors of MSP- 
MS systems. Grimm, Gilfanov & Sunyaev (2002) found that 
Galactic field low-mass X-ray binaries have a scale height 
of ~ 0.410 kpc. Interestingly enough, as shown in Table |3J 
our models over estimate this by almost a factor of two. 
Such an outcome may be another implication that AIC NSs 
receive less momentum at birth than standard NSs formed 
from core-collapse SNe. However, it is not clear that MSPs 
that result from AIC NSs can be linked to an observable 
low-mass X-ray binaries phase (Hurley, Ferrario, Wickra- 
masinghe. Tout & Kiel, in prep). 



5.2 Electron capture supernovae 

Another evolutionary scenario related to NS formation and 
velocity kicks is core-collapse electron capture SNe. This was 
discussed and modelled in Paper I and has also been ac- 
counted for in other population synthesis works (e.g. Ivanova 
et al. 2008). Briefiy, core collapse electron capture SNe are 
thought to arise when electrons are captured onto Mg atoms, 
depleting the electron force in an O-Ne-Mg stellar core of 
sufficient mass (1.4 — 2.5 Mq: Nomoto 1984) which is pro- 
duced by initial progenitor masses in the range of 8 — 12 Mq 
(although this mass range is model dependent: Podsiad- 
lowski et al. 2004). The likelihood that a star born within 



Pulsar Galactic dynamics 21 



the 8 — 12 Mq limit will evolve to have an O-Ne-Mg core 
mass between 1.4 — 2.5 Mq increases if the progenitor is 
able to interact with a companion and lose its outer hydro- 
gen envelope, rather than evolve in an isolated environment 
(Podsaidlowski et al. 2004). Therefore, binary population 
synthesis is ideal for examining the likelihood and outcomes 
of such events. The resultant electron capture SN energy 
yield is low, sufficient to cause the explosion but not enough 
to impart any large velocity to the proto-NS (Kitaura, Janka 
& Hillebrandt 2006). Paper I found that the final MSP spin 
period and spin period derivative parameter space was al- 
tered when electron capture SNe were included. Less pul- 
sar binary systems were disrupted, owing to the small mo- 
mentum imparted during the SN, causing more MSPs to 
be produced. It is also reasonable to expect that including 
electron capture SNe in the binkin models, with SN kicks 
drawn from a distinct distribution with a smaller velocity 
dispersion than for standard NSs, will lead to a reduction in 
the pulsar scale heights. This is a feature that will be fully 
explored in future models so that the impact of the electron 
capture SNe process on binary evolution outcomes and the 
resultant Galactic kinematics of pulsar populations can be 
quantified. 

5.3 MSP-BH systems 

In our models we have assumed no SN velocity kick is given 
to BHs and have found that MSP-BH systems reside close to 
the Galactic plane. If BHs were instead to receive a SN kick 
selected from the same distribution as NSs then it is clear 
that the scale heights of populations including BHs would 
increase (Voss & Tauris 2003). However, it is not so obvi- 
ous that the scale heights would be similar to that of the 
equivalent NS populations (Pfahl et al. 2005). In particular, 
MSP-BH systems (and their progenitors) will be heavier on 
average than MSP-NS systems (and their progenitors) and 
the more massive systems will require a greater momentum 
to reach the same velocities as less massive systems. As such 
MSP-BHs for example, could still have a significant differ- 
ence in their resultant scale height to that of MSP-NSs even 
when both populations receive kicks from the same distribu- 
tion. We would also expect the number of BH binary systems 
to decrease. Most likely it would be the MSP-BH systems 
with a large orbital periods prior to BH formation (systems 
with initial orbital periods greater than ~ 100 days) which 
would be depleted. It is these systems that are not produced 
in the models of Pfahl et al. (2005) who assume SN kicks 
occur on nascent BHs. However, we must bear in mind that 
the final BH masses are calculated assuming that material 
ejected in the SN falls back on to the BH. There is less 
mass-loss associated with BH formation than for NSs and 
this means supernova induced binary disruption is less likely 
during BH formation (in the case of equivalent kick veloci- 
ties). 

We note that when discussing the MSP-BH population 
(or any of our model MSPs) we are defining a rapidly ro- 
tating NS to be an MSP based solely on its spin period. If 
instead we also include consideration of the magnetic field 
strength of these NSs then the nomenclature may be mis- 
leading, especially if we are interested in observable MSPs. 
It turns out that all of the NSs in our model MSP-BH sys- 
tems have magnetic fields residing on, or very close to, the 



assumed bottom mag netic field Umit of 6 x 10^ G (Paper I; 
Zhang & Kojima 2006). Previously (Paper I, Figure |9] and 
Figure [T2|) ■ we have assumed that any NS with a magnetic 
field less than this limit cannot accelerate the electrons in 
its atmosphere to produce the observed pulsar beam and as 
such is not observable as a pulsar. Therefore, if our assump- 
tions regarding accretion on to NSs and how this translates 
to magnetic field decay are correct then we have a lot of trou- 
ble producing observable MSP-BHs. Future observations of 
such systems will help greatly in constraining our evolution- 
ary assumptions. 

5.4 Initial distributions 

In our models we have assumed a maximum birth height 
off the plane, |2maxi|, of either 75 or 150 pc. Consistent with 
Paczynski (1990) and Sun & Han (2004) no significant varia- 
tions of the MSP population scale heights were found when 
varying this parameter. This suggests that the results are 
robust to changes in l^maxil as long as a sensible choice is 
made. The majority of OB star formation has been shown 
by de Wit, Testi, Palla & Zinnecker (2005) to occur within 
\z\ ~ 200 pc of the Galactic plane so choices within this 
range, such as for our models, would seem reasonable. In the 
future it will be interesting to probe the effects of assuming 
a radial dependence in |z;maxi| on the final pulsar popula- 
tion distributions. This may even be tied in with examin- 
ing the effect of assuming bursts of star formation through- 
out the age of the Galaxy and accounting for Galactic arms 
when initiating the birth positions. This final point has pre- 
viously been suggested as an important feature to incorpo- 
rate into population synthesis models (Faucher-Giguere & 
Kaspi, 2007). 

We found that the Yusifov & Kucuk (2004) initial radial 
pulsar birth distribution gave the best fit of our models to 
observations. This distribution was based on observations of 
HII regions within the Galaxy. However, it failed to repro- 
duce the peak of the observable radial distribution - which 
the Paczynski (1990) initial radial distribution succeeded in 
reproducing. Yusifov & Kucuk (2004) recognised that their 
relation is approximate and suggested that a detailed analy- 
sis between models and observations of pulsar velocities and 
Population I stellar positions was required to develop a more 
realistic distribution. We are in a position to do this and as a 
result can suggest that the initial pulsar birth distribution of 
Yusifov & Kucuk (2004) perhaps be shifted towards smaller 
Galactic radii to peak at the inner HII peak (~ 4.0 — 4.5 kpc) 
depicted in Figure 3 of Paladini, Davis & DeZotti (2004). 

5.5 Galactic model potentials 

Even though our favoured model (Model C) utilises the 
Pac90 form of the Galactic gravitational potential we are 
in no way able to distinguish between this and the KG89 
model as a more suitable representation of the Galactic po- 
tential. Both give similar pulsar population scale height re- 
sults which is not surprising given their similarities as shown 
in Figure |4l The form of the XueOS potential is clearly dis- 
tinct from the other two models, especially within the in- 
ner 1 kpc of the Galaxy (where XueOS employ an extrap- 
olation of their measurements), and leads to markedly in- 
creased scale heights. On this basis we do not favour use of 
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the XueOS potential. However, we are not currently in a posi- 
tion to make strong conclusions in this area, especially when 
many previous pulsar, NS and X-ray binary population syn- 
thesis works (such as Paczynski 1990, Lorimer et al. 1993, 
Belczynski, Bulik & Rudak 2002, Sun & Han 2004 and Zuo, 
Li & Liu 2008) have used different Galactic gravitational 
potentials and their results compare well to observations. 
We note that Sun & Han (2004) comment that it is unclear 
whether the Milky Way has a peak in the gravitational po- 
tential at small Galactic radii (as present in the Pac90 and 
KG89 models). 

A possibility in the future is to extend the Galactic 
gravitational potential analysis to consider Modified Newto- 
nian Dynamics (MoND). Such an approach has already been 
taken by Wu et al. (2008), who compare MoND with cold 
dark matter models, and Zuo, Li & Liu (2008) who make 
use of MoND potentials to conduct population synthesis of 
X-ray binaries. 



5.6 Close double compact systems and gamma 
ray bursts 

In this work we have focussed on pulsars and looked in detail 
at MSP systems. However, the models can also be extended 
to explore the formation of close double compact systems 
(NS-NS, BH-BH and NS-BH systems) in detail. The kine- 
matics of these systems is of interest because of their link 
to gamma-ray bursts and, in particular, recent observations 
of the distances at which gamma-ray bursts appear to oc- 
cur from their (assumed) host galaxy (Bloom, Kulkarni & 
Djorgovski 2002). Our combined binpop and binkin mod- 
ules can provide model estimates for the projected distances 
from their host galaxy at which double compact systems co- 
alesce and document the kinematic evolution of these sys- 
tems in general. This will be the focus of a companion paper 
(Kiel, Hurley & Bailes, submitted). 



6 SUMMARY 

We have examined in depth the Galactic dynamics and pop- 
ulation characteristics (owing to stellar, binary and kine- 
matic evolution) of pulsars. Our main findings, reconfirming 
and updating many areas of pulsar evolutionary physics, can 
be summarised as follows (noting that overlap with previous 
work is detailed in Section |4ll: 

• When using a peaked radial distribution for the birth 
locations of binaries, the population of pulsars that arises 
from these binaries also follows a peaked distribution where 
the location of the peak moves inwards in radius by as much 
as 0.5 kpc as the population evolves. Also, compared to the 
birth distribution, the initial shape is preserved inward of 
the peak but the distribution becomes more extended in the 
outer regions. 

• Starting with a uniform initial distribution of binaries 
cannot produce a final pulsar distribution that is peaked 
away from the centre of the Galaxy and therefore does not 
compare well to observations of pulsar locations which indi- 
cate a deficit of pulsars towards the Galactic centre. 

• The form of the Galactic potential does not produce 



significant differences in the final radial distribution of pul- 
sars but can lead to noticeable differences in the calculated 
scale heights of pulsars. 

• As the pulsar population ages the peak of its veloc- 
ity distribution moves to lower velocities. The velocity dis- 
persion of this distribution (assuming a Maxwellian) almost 
halves over a period of 10 Gyr. The shape of the velocity 
distribution is significantly affected by the the inclusion of 
binary evolution - this produces a more sharply peaked dis- 
tribution. 

• Similar to observations we find that the majority of 
standard pulsars are isolated and that these dominate the 
statistics of the pulsar scale height calculations. 

• Isolated pulsars have a greater scale height than binary 
pulsars except in cases where large velocity kicks are ap- 
plied to the population resulting in many isolated pulsars 
being lost from the Galaxy and hence from the scale height 
calculations. 

• Isolated MSPs have greater scale heights than binary 
MSPs (by as much as a factor of two) however, limiting the 
region of the Galaxy considered (in terms of height off the 
plane) does reduce the difference in these scale heights and 
brings them more in line with what observations suggest. 

• We find that 99% of MSPs are in binary systems when 
we only consider SN disruption as a pathway for creating 
isolated MSPs. This does not agree with the observed MSP 
population. If we include a simple ablation model we find 
instead that 64% of MSPs are in binaries which adequately 
matches the observed mix. Furthermore, accounting for ab- 
lation gives similar scale heights for isolated and binary 
MSPs. 

• MSP systems with NS companions can receive large re- 
coil velocities. There is a large scatter in the resulting pecu- 
liar motions of MSP-NS binaries and it is possible for such 
systems to found with low peculiar motion. 

• The scale heights of the MSP-MS and MSP- WD binary 
populations are very similar and follow similar radial dis- 
tributions. These scale heights are larger than that of the 
observed low-mass X-ray binary population in the Galaxy 
(often thought to be the precursors of MSP-MS binaries). 
However, many of the model MSPs in binary systems are 
formed from the accretion-induced collapse of a WD which, 
if given smaller kicks than for standard NSs at birth, would 
reduce the model scale heights. 

• MSPs with WD companions are the most common of 
the binary MSPs. This is followed by MSP-MS, MSP-BH 
and MSP-NS binaries, respectively. 

• Restricting the model MSP population to only include 
MSPs with magnetic fields greater than 6 x lO'^ Gauss dras- 
tically reduces the number of systems and changes the way 
that the population is distributed. This suggests that the un- 
derlying pulsar distribution of the Galaxy may differ greatly 
from the observed sample. 

One future goal of pulsar astronomy is the detection 
of a pulsar orbiting a black hole, and in terms of plac- 
ing constraints upon general relativity a millisecond pul- 
sar in a close orbit around a black hole would be an es- 
pecially exciting observation. We find three distinct evolu- 
tionary pathways which result in the formation of MSP-BH 
systems. These pathways produce three distinct MSP-BH 
populations in terms of orbital period: those with periods of 
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10 day or less, those with periods of about 1 000 days, and 
those with periods of 10 000 day or greater. The short and 
long period populations are the most numerous and only 
the short-period systems are found further than 1 kpc from 
the Galactic plane. We find that owing to the amount of 
material accreted by the MSPs in our model MSP-BH bi- 
naries that the magnetic field decays below 6 x lO'^ Gauss. 
This possibly suggests that we are overestimating the rate 
of accretion-induced magnetic field decay in our evolution 
model - the observation of a MSP-BH binary would confirm 
this possibility. 

We emphasise to the reader that we are not presenting 
any of the models in this paper as a definitive representa- 
tion of the true Galactic pulsar population. The uncertainty 
involved in the many parameters contained within binpop 
and BINKIN does not allow this. Moreover, because we do 
not consider selection effects in our model Galaxy we can- 
not at this stage make definitive comparisons to observations 
as the possibility exists that the observed population may 
be biased in some manner. Lommen et al. (2007) suggest 
that observations of MSPs may preferentially detect binary 
MSPs because the isolated MSPs may be less luminous than 
their binary cousins. The intrinsic luminosity of pulsars is 
not something examined in this body of work. However, it 
will be discussed in detail in future work where selection 
effects are calculated within our upcoming binsfx module 
(Kiel, Bailes & Hurley, in prep). Supplementing our current 
pulsar population synthesis with selection effects will allow 
additional evaluation of the evolutionary codes and their sci- 
entific outcomes. It will also allow us to guide further surveys 
by selecting regions of the sky best suited for the specific pul- 
sar survey and/or telescope of interest. Therefore binsfx will 
provide a powerful tool with which to constrain the theory 
and modelling of stellar, binary and Galactic kinematic evo- 
lution. Further constraints could be placed on binary evolu- 
tion if population synthesis studies are extended to include 
additianal stellar populations and their appropriate selec- 
tion effects. For now, however, we are well on our way to 
producing a comprehensive treatment of pulsar population 
physics. 
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APPENDIX A: 

The time taken for a star to be ablated by a MSP can be 
approximated by taking the irradiated luminosity onto the 
companion, which is of order 

L ~ — ~ 4.4 X 10^^ ergs" \ (Al) 

(Tavani, 1992) and equating AS with the change in binding 

energy of the companion, 

l^bind = —5 , (A2) 



Pulsar Galactic dynamics 

Solving At for the threshold mass of a 0.02 Mq star (Tavani 
1992) gives At = 5.5 Myr, while for a 0.01 Mq star At = 
2.7 Myr. 

This paper has been typeset from a T^JC/ M^^X file prepared 
by the author. 



